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A tutorial discussion of the propagation of waves in random media is presented. In first ap- 
proximation the transport of the multiple scattered waves is given by diffusion theory, but im- 
portant corrections are present. These corrections are calculated with the radiative transfer or 
Schwarzschild-Milne equation, which describes intensity transport at the "mesoscopic" level and 
is derived from the "microscopic" wave equation. A precise treatment of the diffuse intensity is 
derived which automatically includes the effects of boundary layers. Effects such as the enhanced 
backscatter cone and imaging of objects in opaque media are also discussed within this framework. 
In the second part the approach is extended to mesoscopic correlations between multiple scattered 
intensities which arise when scattering is strong. These correlations arise from the underlying wave 
character. The derivation of correlation functions and intensity distribution functions is given and 
experimental data are discussed. Although the focus is on light scattering, the theory is also 
applicable to micro waves, sound waves and non-interacting electrons. 
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I. INTRODUCTION 

Transport of waves through opaque media is a subject of interest in daily hfe. Examples are light transported 
through fog, clouds, milky liquids, white paint, paper, and porcelain, but also electromagnetic waves transported 
through stellar atmospheres and interstellar clouds. Nowadays transport of visible light through human tissue is 
being used as a noninvasive technique to detect, for instance, breast cancer. 

Some basic properties of diffuse light are well known. On a cloud-less day we see the immediate radiation from 
the sun. When a cloud passes in front of the sun, the light first becomes weaker and diffuse. When the cloud has 
become thick enough, the sun becomes invisible; this happens when the cloud thickness is of the order of a mean free 
path. On a cloudy day there is no direct view of the sun, but there is still light: it is diffuse light coming from many 
directions. It has propagated diffusely through the cloud and leaves it in random directions, partly in the direction 
we are looking. 

The study of diffuse wave transport was started by astrophysicists. They wanted to understand how radiation 
created in the center of stars is affected when it traverses an interstellar cloud. Well-known books in this field were 
written by Chandrasekhar (1960) and van de Hulst (1980). For more mundane applications, such as the detection of 
a school of fish using acoustic waves, see Ishimaru (1978). 

It is the purpose of this review to present a comprehensive, self-contained text intended for laboratory applications 
of diffuse wave transport. We explain how the transfer equation follows from the wave equation and how radiative 
corrections can be calculated. Though the approach can be applied to any geometry, we shall focus mainly on slab 
geometries. 

In Sec. H we discuss some general aspects of wave scattering, such as diffusion and Anderson localization. In 



Sec. [II we recall the basic concept of the radiative transfer equation. Next we begin a detailed analysis in section 
[V with the underlying wave equation in the scalar approximation. We explain the notions of t matrix and cross 
section. In Sec. ^ we consider the amplitude or dressed Green's function. We discuss how extinction is related to 
the self-energy. Next we consider propagation of intensity via the Bethe-Salpeter equation. We use the simplest case, 
the ladder approximation, to derive a transport equation equivalent to the radiative transfer equation. This equation 
then describes not only diffusion in the bulk, but also the precise behavior at the boundaries. Various experimentally 
relevant situations are considered in detail: transport in a bulk medium, through very thick (semi-infinite), and through 
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finite slabs; the enhanced backscattcr cone; exact solutions of the Schwarzschild-Milne equation; the regime of large 
index-mismatch; semiballistic transport; and imaging of objects hidden in opaque media (sections VI-XIII). 

With these ingredients we consider in the second part correlations of intensities. This leads to three kinds of 
observable correlation fimctions. the angular correlation function Ci, the correlation of total transmission C2, and the 
correlation of the conductance C3. We explain how these functions are calculated and next we consider correlations 
between three intensities (third cumulants) and the full distribution of intensities. 

A. Length scales 

The diffuse regime is formulated in the three inequalities (Kaveh, 1991), 



where A is the wavelength. C. the mean free path, L the sample size. Labs the absorption length (for optical systems), 
and Line the incoherence length (for electronic systems). The first inequality ensures that localization effects (see 
below) are small; the second inequality implies that many scatterings occur if the wave transverses the system; the 
third inequality ensures that not all radiation is absorbed. 

The description of radiation transport can occur on roughly three length scales: 

-Macroscopic: On scales much larger than the mean free path the average intensity satisfies a diffusion equation. 
The diffusion coefficient D enters as a system parameter that has to be calculated on mesoscopic length scales. 

-Mesoscopic: On length scales of the mean free path ^, the problem is described by the radiative transfer equation 
or Schwarzschild-Milne equation. This is the Boltzmann-type equation for the system. At this level one needs as 
input the mean free path ^ and the speed of transport v, which should be derived from microscopies. In the diffusive 
regime this approach leads to the diffusion coefficient D = v£/3. 

-Microscopic: The appropriate wave equation, such as Maxwell's equations, the Schrodinger equation, or an acoustic 
wave equation, is used on this length scale. The precise locations and shapes of scatterers are assumed to be known. 
Together with the wave nature they determine the interference effects of scattered waves. In light-scattering systems 
the scatterers often have a size in the micron regime, comparable to the wavelength A, which could lead to important 
resonance effects. The mesoscopic or Boltzmann description follows by considering the so-called ladder diagrams. 

The microscopic approach will be the starting point for this review. Fundamental quantities, such as the self-energy 
and the Hikami box, can only be calculated on this level. The drawback of the microscopic approach is that it is 
too detailed. In practice the precise shape and positions of the scatterers often are not known and a mesoscopic or 
macroscopic description is necessary. 

B. Weak localization, closed paths and the backscatter cone. 

The diffusion equation is a classical equation that fully neglects interference effects inherent to wave propagation. 
At this level of description there is no difference between diffusion of particles and of wave intensity. Whereas 
a transmission pattern of monochromatic (laser) light through an opaque medium is known to consist of speckles 
(bright spots on a dark background), the diffusion equation and the radiative transfer equation only describe the 
average intensity. 

The wave nature of light immediately leads to a reduction in transmission due to interference effects. Following 
Rayleigh (1880) we suppose that the transmission amplitude E for a certain experiment is composed of many terms 
Cp arising from physically different interference paths p. 

Some paths have closed loops, i.e. loops that return to the same scatterer. When such loops contain two or more 
common intermediate scatterers, they can be traversed in two directions. Let us consider one of such loops Cp, and 
denote the contribution of the second, reversed loop by Dp. Summing over all paths we have 



A<^<L<Labs,ii: 



inc t 



(1) 




(2) 



p 



The intensity is 



= Y,ic;cp + d;Dp) + Y.{c;Dp + D;Cp) + ^ {c; + D;){Cp, + Dp,). 



(3) 
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When applied to electrons, quantum mechanics tells us that the C*C and D*D terms are probabilities, while C*D 
and D*C terms are interference contributions. Naively, one expects the second and third term to be small. Thus in 
Boltzmann theory only probabilities are taken into account, that is to say, the first term. 

However, if there is time-reversal invariance, the second term X^p(C'*Dp + D*Cp) will be equally large: there is a 
factor of 2 for each closed loop. As a result, for the wave intensity there is a larger probability of return. 

In optics there is a simply measurable effect due to this, namely the enhanced backscatter cone. When the 
incoming and outgoing light are in exactly the same direction, the light path may be considered closed. As predicted 
by Barabanenkov (1973) and detected by several groups (Kuga and Ishimaru, 1985; van Albada and Lagendijk, 
1985; Wolf and Maret, 1985) the average intensity in the backscatter direction has a small cone of height almost 
one and angular width S6 X/i. This observation has given an enormous push to the research of weak localization 
phenomena in optics. 

This effect is a so-called weak localization effect. These effects occur if X/£ <^ 1 and are precursor of the so-called 
strong localization effects which occur ii X/£ ^ 1. These effects are also known as "mesoscopic," indicating length 
scales between macroscopic (the diffusion equation) and microscopic (individual scattering events). Indeed, in electron 
systems these effects only show up in rather small samples due to inelastic scattering. For optical systems there is no 
such size restriction. The study of mesoscopic effects started in the field of electronic systems, and we shall use many 
results derived there. 

For electronic systems weak localization effects were first analyzed by Altshuler, Aronov and Spivak (1981) in order 
to explain the Sharvin-Sharvin fluctuations of a resistance as function of the applied magnetic field (Sharvin and 
Sharvin, 1981). These so-called "magneto-fingerprints" show a seemingly chaotic conductance as the applied field is 
varied, but are perfectly reproducible as they are solely determined by the location of the scatterers. (Only at very 
low temperatures electron scattering is dominated by impurity scattering; at higher temperature scattering mainly 
arises from phonons.) 

C. Anderson localization 

If scattering is very strong, the "weak localization" effects become of order unity, i.e., ^ ~ A ( A/£ is of the order 
of 0.001-0.01 for visible light in standard laboratory situations). According to the criterion of loffe and Regel (1960) 
the diffusion constant will tend to zero at that point. The scattering process that causes the cinhancxid backscatter 
also reduces the diffusion constant. The amplitudes that make up the intensity, split and one of the amplitudes visits 
some of the scatters in reverse sequence, thus forming a loop. The diffusion constant is lowered by these processes. 
(The loop is somewhat similar to the renormalization of the effective mass in field theories.) If scattering is increased, 
the contribution of these loops becomes more and more important. It is not simply that the diffusion constant tends 
linearly to zero as scattering increases. Rather, the return probability of the intensity becomes higher and higher, 
reducing the diffusion constant in a stronger fashion. The diffusion constant can thus become zero at finite scatterer 
strength, meaning that the wave can no longer escape from its original region in space. This is the transition to the 
well-known Anderson localization (Anderson, 1958). 

In the Anderson localization regime there are only localized states. In the dclocalizcd regime there exist extended 
states, responsible for diffusion. In the localized regime the intensity typically decays exponentially over one local- 
ization length, whereas in the diffuse regime the wave function extends up to infinity. There are thus two different 
phases, the diffuse, metallic regime and the localized, insulating regime. In three dimensions a phase transition from 
the extended to the localized state can occur. In one and two dimensions the states are always localized, provided 
the waves are non-interacting (for electrons also spin-orbit scattering has to be absent). Yet for a finite sample the 
localization length can be much larger than the system size, in which case the states appear to be extended and the 
conductance does not vanish. Note that the localization is solely the result of the interference of the waves scattered by 
the disorder. (This is not the only scenario for a metal-insulator transition in electron systems. Due to their fermionic 
nature and interactions, electron systems allow for a whole range of possible transitions between the conducting and 
insulating regimes; see Mott (1990)). Anderson localization is also called strong localization. 

The precise behavior near and at the transition is not fully understood. The standard diagrammatic perturbation 
technique, which we use, works well for the description of diffusion and low order corrections, but it is not suited for 
the study of the transition. Therefore various techniques have been developed to study behavior near the transition 
and the phase transition itself. 

An important step was the scaling theory of localization put forward by Abrahams, Anderson, Licciardello and 
Ramakrishnan (1979), which states that near the localization the only parameter of importance is the dimensionless 
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conductance g (the conductance measured in units of e'^/h). The scaUng of g as a function of sample size was 
studied earlier (Thouless, 1974; Thouless, 1977; Wegner, 1976). Abrahams et al. extended those ideas and derived 
renormalization-group equations. They concluded that in one and two dimensions there is no real phase transition; 
the states are always localized. In three dimensions a phase transition can occur. 

Classical diffusive transport is described by the ladder diagrams. The so-called maximally crossed diagrams are 
the next most important diagrams in the diffuse regime. They describe the leading interference terms responsible 
for the backscattcr cone and reduction of the diffusion constant. VoUhardt and Wolfle summed sclf-consistcntly 
combinations of the ladder diagrams and the maximally crossed ones in the diffuse regime (VoUhardt and Wolfle, 
1980a; VoUhardt and Wolfle, 19806; VoUhardt and Wolfle, 1992). They found that the diffusion constant vanishes at 
strong enough scatterer strength, thus providing a microscopic picture of the Anderson transition. As the maximally 
crossed diagrams yield loops of intensity, the approach can be seen as a self-consistent one-loop summation. Although 
a self-consistent approach will certainly not include all diagrams of higher order, the method works flne even close the 
transition (Kroha, 1990; Kroha, Kopp and Wolfe, 1990), because the first few higher order correction terms vanish. 

Another approach is to perform an exact average over the disorder within a field-theoretic approach. Next, one 
integrates over the fast fluctuations and is left with the slow variables of the system. This technique yields the so-called 
nonlinear sigma model (Wegner, 1979; Hikami, 1981; Efetov, 1983; Altshuler, Kravtsov and Lcrncr, 1991). With the 
resulting action it is possible to generate systematically all corrections to the diffusion process (in 2 -|- e dimensions) , 
allowing for an explicit foundation of the scaling theory. 

A recent approach is that of random matrix theory. The basic assumption here is that the total scattering matrix of 
the system, although very complicated, can be described just by random matrix elements respecting the symmetries 
of the system. Surprisingly, this method works well, and just from the assumption of the random ensemble it predicts 
many features of the systems correctly and in a simple way. Its is, unfortunately, only applicable to quasi one- 
dimensional situations. An overview of the theory and applications of random matrix theory is given by Stone, Mello, 
MuttaUb and Pichard (1991) and Beenakker (1992). 

Although some years ago there was hope that the Anderson transition might soon be reached for light scattered 
on disordered samples, it was not observed when this review was written. One did observe that in time resolved 
transmission measurements the average transmission time became very long. This was interpreted as an indication 
that the mean free path was very small, which would mean that the light was close to localization. It turned out, 
however, that the long transmission times arise since the light spends much time inside the scatterers. This meant 
that the Anderson localization was still out of reach (van Albada, van Tiggelen, Lagendijk and Tip, 1991). 

D. Correlation of different diffusons 

Another interaction effect is the interference of a diffuse intensity with another that has, for instance, a different 
frequency or position. In this work wc shall concentrate on such processes which, lead to correlations in, for instance, 
the transmitted beams. There are advantages to studying the correlations above the loop-effects. The correlations can 
be measured more accurately and easily in experiments than can renormalization effects. Secondly, it is an interesting 
feature of optical systems that there are three different transmission measurements: 

- Coming in with a monochromatic beam in one direction (this we call "channel a" ) one can measure the intensity 

in the outgoing direction b and define the angular transmission coejjicient Tat- Its correlation function is called the 
Ci correlation and is of order unity, describing the large intensity difference between dark and bright transmission 
spots. 

- One can also measure all outgoing light. In practice one uses an integrating sphere. This leads to the measurement 
of the total transmission, 



Its correlation function is called the C2 correlation function. 

- Finally, one can also add the results of coming in from all possible directions, either by repeating the experiment 
under many different incoming angles or by using diffuse incoming light. This leads to a quantity 




(4) 



b 




(5) 



a 



ab 
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In analogy with electronic systems it is called the conductance. Its fluctuations are called C3 fluctuations in optics. 
In electronics (where, for instance, the magnetic field is varied), they are called universal conductance fluctuations 
(UCF). Notice that these fluctuations are not temporal but static, since the scatterers are fixed. Both C2 and C3 
correlations are low order corrections in 1/g, which acts as the small parameter giving the relative strength of the 
interference effects. 

There are various books and proceedings on localization and mcsoscopics in electron systems (Nagaoka and 
Fukuyama, 1982; Ando and Fukuyama, 1988; Kramer and Schon, 1990; van Haeringen and Lenstra, 1990; van Haerin- 
gen and Lenstra, 1991; Altshuler, Lee and Webb, 1991; Hanke and Kopaev, 1992) and classical waves (Kirkpatrick, 
1985; Sheng, 1990; Soukouhs, 1993). Many aspects of the Anderson transition have been discussed by Lee and 
Ramakrishnan (1985). 



II. MACROSCOPICS: THE DIFFUSION APPROXIMATION 

A. Transmission through a slab and Ohm's law 

Consider the propagation of light through a slab of thickness L ("plane parallel geometry"). As indicated in Fig. |l| 
a plane wave impinges on the surface of an opaque medium at an angle 9a ■ The index of refraction of the medium is 
no, and that of the surroundings is ni. The ratio of the indices of refraction, 

m= — , (6) 

ni 

is larger than unity if a dry substance is placed in air (ni w 1). If the substance is placed in a liquid, one may have 
TO < 1. If TO 7^ 1 the refracted beam inside the medium will have an angle 9'^ different from 9a. 
Due to multiple scattering the incoming beam decays exponentially, 

/™(z) = /oe-^/^-^''", (7) 

where the cosine is a geometrical factor expressing the total path length in terms of z. This exponential decay is 
known as Lambert-Beer law. It describes the decay of the direct sunlight in clouds or of headlights of cars in fog. 
The light source becomes invisible when the thickness of the diffuse medium is greater than one mean free path. The 
decay of direct light intensity occurs because it is transformed into diffuse light. 

The scattering of the incoming beam into diffuse light occurs in a skin layer with characteristic thickness of one 
mean free path. Later on this will be discussed in greater length when we solve the Schwarzschild-Milne equation. 
In the diffusion approach one simply assumes that the incoming beam is partly refiected and partly converted into 
a diffuse beam in the skin layer. Next one assumes that effectively diffuse intensity enters the system in a trapping 
plane located at distance zq ^ i outside the scattering medium. The value of zq is phenomenological, but a precise 
analysis of the Schwarzschild-Milne equation reveals that this picture is valid. One usually takes zq — 0.7104£, the 
exact value for isotropic point scattering of scalar waves, van de Hulst (1980) investigated many possible cross sections 
and numerically always found values very close to this. Amic, Luck and Nieuwenhuizen (1996) considered the limit 
of very strong forward scattering and found that zq lies only 1.1% above the quoted value for isotropic scattering. 

Once the light has entered the bulk, the diffuse intensity obeys the diffusion equation 

^^I{r,t)^DV'l{r,t). (8) 

In the steady state the time-derivative vanishes and the diffusion coefficient plays no role. For a slab geometry with 
plane waves there is no x,y dependence, and we have to solve I"{z) = 0. The trapping plane boundary conditions 
are /(— 20) = loi I{L + zo) = 0. The solution is a linear function of z. 

The transmitted intensity is essentially equal to this expression at z = L (at a distance zq ^ ^ from the trapping plane 
at L + Zq) and equal to Zq times the derivative at that point; this estimate will turn out to be qualitatively correct 

I{Z = L) _ Zq dl{z) 
J — — : \^=L+zo- 



lo lo dz 

Zq Zq 



L + 2zo L 



(10) 
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Without having paid attention to details we find the behavior T ^ IjL^ equivalent to Ohm's law for conductors. 
Generally a conductor has conductance jh per channel. Here there are Aj}? channels, of which a fraction T is 
open. This yields the conductance 

^ A e^lA , , 



B. Diffusion propagator for slabs 

In the presence of absorption the diffusion equation reads 

9f/(r, i) = D\7^I(r, t) - DK^I{r, t) + ^(r, t), (12) 

where is a source term, the second term on the right-hand side describes absorption, and the absorption length is 
Labs = In a bulk system the solution with initial conditions /(r, 0) = 5{v) and 5 = reads 

/(r,i) = (47rDi)"^^^e-3TC--DK't. (13) 

Its Fourier-Laplace transform 

/(q, j dVe-^'i '' die-*"*/(r, t), (14) 

takes the form 

For HI — these expressions diverge in the limit q, O ^ 0. This divergency is called the "diffusion pole" or "diffuson" . 

For a slab geometry there will only be translational invariance in the p — {x,y) plane. Suppose we have a source 
S{r) = S{x,y)S{z — z'). Denoting the perpendicular wave vector by q_L we have for diffusion from z' to z 

I"{z, z') = M'l{z, z') - ^^Siz - z'), (16) 



with the "mass" M defined by 



M^ = q^^ + n^ + ^^. (17) 



M is the inverse depth at which a given intensity contribution ■i/'aV^fe of an incoming beam -ip — has decayed by 

a factor 1/e, due to spatial de-phasing of the amplitudes (encoded in q^), temporal de-phasing (expressed by SI) or 
absorption (expressed by k.) 

By realizing that the diffusion equation Eq. (^) is just a wave equation with complex frequency, one sees that the 
solutions to the diffusion equation are a linear combination of hyperbolic sines and cosines. (The solution can also be 
obtained using the method of "image charges" known from electrostatics.) The diffuse intensity propagator has the 
form (Zhu, Pine and Weitz, 1991; Lisyansky and Livdan, 1992) 

I(z z'- -n)-^^^-^^ [sinhMz< MzoCOshMz<][sinhM(L-z>) + MzoCoshA/(i-z >)] 
(z,z,(l±, ) - — ^ (M + M^zl) sinh ML + 2M^zo cosh ML ' ^ ^ 

where 

z< = min(z, z'), z> = max(z, z'). (19) 



In the stationary limit in the absence of absorption (fi = /t = 0, q^ = ^ A/ = 0), Eq. (21) reduces to a tent-shaped 
function. 
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(20) 



The propagator describes the diffuse propagation from one point in the slab to another. With roughly equal indices 
of refraction inside and outside the sample, the extrapolation length zq is a few mean free paths and thus the terms 
involving zq yield contributions of order £/L. For optically thick samples {L £), this is negligible and one has 

Tf ' OA ^(q^) sinh(Mz<)sinh(Af£-Afz>) 

I{z, z;ci^;n) = — Msinh(AfL) ' ^^'^ 

The diffusion equation does not hold if the intensity gradient is steep, i.e., if qi ^ I. It can therefore not describe 
properly the diffuse intensity near the surface of the medium. Nevertheless, in a heuristic manner one often makes the 
following assumptions in the diffusion approach (Ishimaru, 1978): (l)The diffuse intensity from an outside plane- wave 
source is assumed to be given by substituting z = £ in Eq. (|l8|), as if all diffuse intensity originated from this z = £ 
plane. Likewise, the coupling outwards is obtained by taking z' — L — £. In this way two extra propagators, the 
incoming (z < 0, < z' < L) and the outgoing (Q < z < L, z' > L) are constructed from the internal diffuson 
(0 < < L). (2) The boundary conditions arc (Zhu et ai, 1991) 

H^, - zo dj{z, z'- M)|,^„+^^_ . (22) 

The form Eq. ( p^ ) fulfills these conditions. The extrapolation length zq is determined by the reflectance at the surface. 
In the following we calculate it precisely. 



III. MESOSCOPICS: THE RADIATIVE TRANSFER EQUATION 

The study of multiple light scattering was initiated in astrophysics with the goal to derive, on the basis of energy 
conservation, the radiative transfer equation. This is the "Boltzmann" transport equation of the problem (i.e., the 
mesoscopic balance equation that neglects all memory effects). It has been solved in particular for slabs (plane-parallel 
geometries). This approach has been described in the books by Chandrasekhar (1960) and van de Hulst (1980). It 
will be shown below that the radiative transfer equation can also be derived from the ladder approximation to the 
Bethe-Salpeter equation. In other words: there exists a microscopic derivation of the radiative transfer equation. 
Once this is shown, more details can be incorporated in the microscopies and more subtle effects, such as backscatter 
and correlations, can be derived microscopically in a way closely related to the radiative transfer equation. 



A. Specific intensity 

The specific intensity X(r, n) = I(r, 6*, (p) is defined as the radiation density emitted at position r in direction 
n = (sin 6* cos (/), sin 6* sin(/), cos6') in a system with density n of scatterers. Let AU be the radiation energy in a given 
frequency interval {lu — ^IS.uj,uj + 5 Aw), transported through a surface dcr in directions lying within a solid angle dn 
centered around n, during a time interval dt. This energy is related to the specific intensity as (Chandrasekhar, 1960) 

d[/ =Xcose Aw do- dn dt, (23) 

where Q is the angle between the director n of the emitted radiation and the normal of dcr. The energy depends in 
general on the position r = {x, y, z), the direction n, the frequency w, and the time t. 

We are mainly interested in stationary, monochromatic situations for a slab with axial symmetry. Then X depends 
only on z and = cos 0, where 9 is the angle between the z axis and the direction of emitted radiation. We consider 
propagation in a medium with density n of rotationally symmetric scatterers with extinction cross section aex and 
phase function p{n,n') — p{n ■ n') = p(cosG). Thus p{n ■ n') dn' dn is the fraction of radiation entering inside a 
narrow cone of width dn' around incoming direction n' and leaving inside a narrow cone dn around n. For spherically 
symmetric scatterers p is a function of n' • n = cos0. 

If radiation propagates over a distance ds, there is a loss of intensity due to scattering into other directions and 
due to absorption (both are incorporated in aex)-, 

dX = —ncTexXds. (24) 
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There is also a gain term naexSA.s with the source function J ^ 



j{v-e,cj^)^ j p{e,cj^-e'<p'f-^^^^x{v-e\ci^'). (25) 

S describes the radiation arriving at r in direction n' and scattered there in direction n. The radiative transfer 
equation expresses the net effect of the gain-loss mechanism 

^f = J-I. (26) 
nUex ds 

Here it is derived from phenomenological considerations, later we will provide a microscopic derivation. 

The loss term —I leads to the Lambert-Beer law. Indeed, for a unidirectional beam X(0, n) = /o(5(n — no), simple 
integration yield I(r, n) = Io5{n — no) exp(— r/£sc)+scattered. Which is again the Lambert-Beer law with scattering 
mean free path 

4c = (27) 



B. Slab geometry 

For homogeneous illumination of a slab < z < L, physical quantities depend only on the depth z. It is useful to 
introduce the "optical depth" 

r-f. (28) 



The optical thickness of the slab is then 



6=-^. (29) 



Let 9 be the angle between the direction of radiation and the positive z axis, and <j) its angle with respect to the 
positive X axis. This allows us to introduce the dimensionlcss form of the radiative transfer equation, 

M^^^5^ = ^(r,M,0)-I(T, /.,</-), (30) 

where 

/i = cos 9. (31) 

For a plane wave with intensity /o incident under an angle xia{9a, 4>a) on the interface z = 0, the boundary condition 
is X(0, ^, 0)= IoS{fi - tJia)&{4> - 4>a), where //q = cos 6*0. > 0. 



1. Isotropic scattering 

The radiative transfer equation (^0|) can be written as an integral equation. For a slab this is usually called the 
Schwarzschild-Milne equation or, for short, the Milne equation. Let us consider a semi-infinite medium. For /i < 
Eq. ( ^ ) yields for the radiation in the —z direction (backward direction) 

r°° dr' 
J(r,M, </-)-/ :7(r',Ai,0)e-(^-^)/l^l— , (32) 

while for /i > the specific intensity in the -\-z direction satisfies 

X(r, /i, (/.) = J(0, /i, ct^y-^'^ + r J{t\ 0)e-(^-^')/^ — . (33) 

Jo M 
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Of special interest is the case of isotropic scattering or s-wave scattering for which 

p(cose) = l. (34) 

Isotropic scattering occurs for electron scattering from small impurities. For isotropic scattering J7 does not depend 
on /i and (j). It is common to introduce the dimensionless intensity 

r(T) = -iyd/x'd0'X(r,^V), (35) 
which here equals AttJ'/Iq . Combining the last two equations yields for a plane wave incident in direction (fia, 4'a) 

r(r) = e-^/^" + r dr' C ^T{t'). (36) 

Jo Jo 2^ 

This Boltzmann equation will also be found from the ladder approximation to the Bethe-Salpeter equation, which 
provides a microscopic foundation. Corrections to the ladder approximation then yield the limit of validity of the 
Boltzmann approach. 

The precise solution of Eq. (|3^) can be obtained numerically. See van de Hulst (1980) and Kagiwada, Kalaba and 
Ueno (1975) for details. We have plotted it in Fig. ^ for a relatively thin slab (L = Ai), using the data from Table 
17 in van de Hulst (1980). The albedo is unity, and no index mismatch. The form of the solution near the incoming 
surface is quite different for the three cases drawn: fj,a = 1, or perpendicular incidence; /ia = 0.1 (angle with the z 
axis is 84 degrees); and diffuse incidence uniformly distributed over all angles. At the outgoing surface all solutions 
are alike: there is only a small deviation from the straight line crossing zero at L + zq. 



2. Anisotropic scattering and Rayleigh scattering 

For acoustic or electromagnetic waves scattered from particles much smaller than the wavelength there is no s-wave 
scattering. Instead one has to leading order the dipole-dipole or p-wave scattering. It is expressed by the phase 
function for Rayleigh scattering, 

p(cose) = ^(1 + cos^e). (37) 

Scattering from spherically symmetric particles is usually anisotropic, but cylindrically symmetric with respect to 
the incoming direction. For arbitrarily shaped scatterers this symmetry holds only after averaging over their possible 
orientations. In such situations the phase function p(cos0) depends only on — 0', and the average over can be 
carried out. This leads to the projected phase function, with /i = cos 9, 

Po(M,/i') = / ^p{0<t>;O'^') = J ^^pie^;0'^'). (38) 

For Rayleigh scattering one has, using the equality cos 9 = n • n' = sin sin 9' cos(0 — (j)') + cos 9 cos 9', 

Po(m, m') = ^(3 - - /i" + 3A*'/i"). (39) 
Because of the form of the integral, it is useful to introduce 

r(r) = ^ df,I{r, Ai) A(r) ^ ^ d/iAi'X(r, f,). (40) 

The radiative transfer equation yields the coupled integral equations 

r(r) = e-^/^" + i^E^ - ^E^) * F + {^E^ - ^E,) * A 

A(r) = M^c-/'- + {^E3 ^E,) * F + {^E, - ^E^) * A. (41) 
Here we introduced the exponential integrals Ek, 



11 



Ja M Ji r 



(42) 



and the product 

{E*f){r)^ I dT'E{\T-r'\)f{T'). (43) 



Eqs. ( p4| ) involves two coupled functions and represent the simplest extension of isotropic scattering. They have been 
analyzed by van de Hulst (1980). 



3. The transport mean free path and the absorption length 

We have discussed how unscattered intensity decays exponentially as a function of the distance from the source. 
This occurs because more and more light is scattered out of the direction of the beam. The characteristic distance 
between two scattering events is the scattering mean free path £sc- Now suppose that scattering is rather ineffective, 
so that at each scattering the direction of radiation is not changed much. The diffusion constant must then be large. 
In other words, the factor £ in the identity D = \vl cannot be the scattering mean free path. Intuitively one expects 
that £ is the distance over which also the direction of radiation gets lost. This length scale is called the transport 
mean free path. We show now how it follows from the radiative transfer equation. 

The time dependent radiative transfer equation has the form 



d f dn' 

isc^2:(r,n,t) +4cn • VX(r,n,i) = / — p(n, n')I(r, n', <) - X(r, n, i). 



(44) 



where tsc is the mean time between two scatterings. We can now introduce the local radiation density / and the local 
current density J as 



I{r,t) — J dnX(r, n, i), J{r,t) = y~ J dnl(r,n,t)n. 
Integration of Eq. (H) yields the continuity equation 



(45) 



ft/(r,i) + V.J(r,t) = -l-^/(r,t). (46) 

^SC 

Here we introduce the albedo (from the Latm albus^ white), the whiteness of the scatterer, 

a = l^p{n,n'). (47) 

For a = 1 there is no absorption. If we multiply Eq. ( ^4[ ) by n and integrate over n we obtain 

^atJ + ^J + 4cV- / dnl(r, n, i)nn = / dndn'p(n, n')X(r, n', i)n. (48) 

If scattering is (on the average) spherically symmetric, the integral / dnp(n, n') n can only be proportional to n'. If 
one takes the inner product with n', one finds the average cosine of the scattering angle, 

/dn P dn d 

-£p{n, n') n • n' = j ^p(cos 6) cos 9 = j ^ ^p(a*)a'- (49) 

Therefore Eq. ( ^ ) can be written as 

t r 

-^StJ + -^(l-(cose))J--4cV- / dnl(r,n,i)nn. (50) 

The right-hand side depends on / and J. If one assumes that the intensity distribution is almost isotropic, then the 
current is much smaller than the density. This allows one to make the approximation (Ishimaru, 1978) 
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3f 

J(r, n, t) « /(r, t) + -^n • J(r, t) + • • • , (51) 
since higher order terms are smaller. Under these approximations it follows that 

f^dtJ + (1 - (cose))^J ^ -%V/. (52) 

^sc ^sc "J 

For processes that change slowly in time we thus find 

J(r,<) = -DV/(r,t). (53) 



When inserted in the continuity equation, Eq. (53) leads to the wanted diffusion equation for the density. Thus under 
the above assumptions the radiative transfer equation leads to the diffusion equation, 

dJ(r,t)=DW^I{r,t)~^I{r,t) 

^SC 

= DV^I{r,t)-DK^I{r,t). (54) 

The diffusion constant D is given by 

D=——^^-——^\vetr, (55) 
3tsc(l - (cose)) 3 

where v = Iscl^sc is the transport speed. We wish to stress that in principle tsc has two contributions: the time to 
travel from one scatterer to the next and the dwell time spent in the neighborhood of one scatterer(van Albada et at, 
1991). In Eq. ( [5^ ) the transport mean free path occurs 

(tr - l^-^isc (56) 

1 — (cos B) 

It is the mean distance after which the direction of radiation gets lost. For strongly forward scattering (cos 6) will 
be close to unity so that the transport mean free path becomes large. In this excercise we also found an explicit 
expression for the absorption length Labs and the inverse absorption length k 

"-^'■'^-n^^^w^y ^^^^ 



C. Injection depth and the improved diffusion approximation 

In the previous section we avoided the problem of how an incoming plane wave becomes diffusive inside the medium. 
We introduced a "trapping plane" at the injection depth zq ~ £. The precise statement is that the solution of the 
transfer equation, or Milne equation, of a semi-infinite medium starting at z = behaves as I{z) = z + zq for z ^ £; 
formally this expression vanishes at — zq outside the medium. This can also be expressed as 

/(O) = zo/'(0). (58) 

If there is an index mismatch between the system and its surroundings, the walls of the system will partially reflect 
the light. Therefore the light will remain longer in the system. The transmission coefficient will become smaller by a 
factor of order unity. This effect is of practical importance, as the scattering medium usually has a different index of 
refraction from its surroundings, usually air or liquid. Only in the special case of index matched liquids the mismatch 
is minimized. 

The first to point out the importance of internal reflections were Lagendijk, Vreeker and de Vries (1989). They also 
noted that zq changes. For a one-dimensional medium they give the expression 

zo - l^e^a, (59) 
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in which R is the mean reflection coefiicient. 

For the three-dimensional situation, an analogous result was worked out by Zhu et al. (1991): For a system in which 
only the z dependence is relevant, Eq. (|5l|), relating the specific intensity X to the radiation density / and the current 
density J, reads 



X(z, ^) = I{z) + -fiJz{z), 



(60) 



with velocity v = Isc/tsc, where tg^ is the scattering time. It follows that 



/I n2-K 



d(j)2{z, ii). 



From this one reads ofi^ that the radiation current per unit of solid angle dQ = d/i dcf) equals 

^ = -a'-:| = ^WW + 3a^V(z)}. 

The total radiation at depth z = in the positive z direction is thus 



(61) 



(62) 



V 

47r 



/id/i27rX(0, /i) 







^/(0) + ij,(0) 



(63) 



where the diffusive current (53) has been inserted. 



In the absence of internal reflections dJ^/dfi must vanish for all /i > 0. If we impose that J^^ vanishes and if we 



compare with Eq. (58), we obtain 



(64) 



For isotropic scattering this expression is not far from the exact value zq — O.llQAM that follows from the radiative 
transfer equation, see Sec. 

Let us now assume that the refractive index of the scattering medium uq diflers from that of its surroundings ni. 
The ratio 



no 
ni 



(65) 



exceeds unity for a dry medium in air, but can be smaller than unity if the medium is in between glass plates with a 
high refractive index. In both cases internal reflections will appear at the interface. The reflection coefficient is given 
by 



R{p) 



ji — \/m ^ — 1 — /i^ 



1-Ai^ 



(66) 



i?(/i) equals unity in the case of total reflection, namely, when the argument of the square root becomes negative. We 
can now calculate the current that is internally reflected at the interface z = 0^: 



3 



rcfl 



d^ / d</>i?(/x) 



dJ, 
dVL 



d/i 

-1^ 



i?(M) (1^^/(0) + 3mV,(0)) 



where 



Ci = d^^R{^i), C2 = / d/i/i" 
Jo Jo 



(67) 



(68) 
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By equating this angular average to the angular average (|6^), Zhu, Pine and Weitz derive Eq. (jS^) with 

2 1 + 3C2 „ 

It was later pointed out by Nieuwenhuizen and Luck (1993) that this result becomes exact in the limit of large index 
mismatch (to ^ or m ^ oo). The reason is that then the interfaces are good mirrors, so that the outer world is not 
seen and hence also close to the mirrors the intensity is diffusive. 

In the derivation of Eq. (|6^) the phase function has only been used in the expression for the diffusion coefficient. 
This led to the occurrence of the transport mean free path. One would therefore expect that Eq. (^9|) remains valid for 
arbitrary anisotropic scattering in the limit of large index mismatch. This can be explained as follows: The interfaces 
act as good mirrors. Therefore many scatterings also occur close to the wall before the radiation can exit the medium. 



After many scatterings the radiation has become isotropic anyhow. Evidence for this point will be given in Sec. X.D 



IV. MICROSCOPICS: WAVE EQUATIONS, T MATRIX, AND CROSS SECTIONS 
A. Schrodinger and scalar wave equations 

Light scattering is described by the Maxwell equations. The vector character of the amplitudes leads to a tensor 
character of the intensity. In regard to sound waves or Schrodinger waves, this introduces extra complications, that 
we do not address here. For discussions of the vector case see, e.g., van de Hulst (1980), MacKintosh and John 
(1988), Peters (1992), and Amic, Luck and Nieuwenhuizen (1997). The vector character is especially important in 
case of multiple scattering of light in a Faraday active medium in the presence of a magnetic field. For fundamental 
descriptions, see MacKintosh and John (1989) and van Tiggelen, Maynard and Nieuwenhuizen (1996). This field 
also includes the so-called photonic Hall effect, a light current that arises, under these conditions, perpendicular 
to the magnetic field. From studies of the diffuse intensity Nieuwenhuizen (1993) imagined such a mechanism to 
exist. This was indeed shown by van Tiggelen (1995), and confirmed experimentally by Rikken and van Tiggelen 
(1996). Another interesting application is multiple scattering of radiation emitted by a relativistic charged particle in 
a random environment, see Gevorkian (1998) for diffusion, and Gevorkian and Nieuwenhuizen (1998) for enhanced 
backscatter. 

Acoustic waves and spinless electrons are described by scalar waves and this turns out to be a good approximation 
for light as well. The Schrodinger equation reads in units in which /2m — 1 

-V^* + = Em. (70) 

As potential we choose a set of point scatterers. 

V{v)^-Y,u5{v-n,). (71) 

i 

Here —u is the bare scattering strength and Ri are the locations of the scatterers. 
Acoustic waves are described by the classical wave equation 

where e is the normalized local mass density and c is the speed of sound in a medium with e = 1 . We shall apply the 
same equation to light, thereby neglecting its vector character, e is then the dielectric constant, and c is the speed 
of light in vacuum where e = 1. For monochromatic waves ^'(r, i) = ^'(r)e'"* the time derivatives are replaced by 
frequency iuj. In practice this applies to a stationary experiment with a monochromatic beam. For classical waves we 
can also introduce point scatterers by setting £(r) = 1 + ai5(r — R^), where a is the polarizability of the scattering 
region that we approximate by a point. The Schrodinger equation and the scalar wave equation then both take the 
form 

- V2^'-uX;.<5(r-Ri)* = /fc^*, (73) 

where / " ~ constant in the Schrodinger equation; k = \/E 

\^ u = ak^ in the classical wave equation; k — uj/c. 

Many static results can be derived without specifying which kind of waves are discussed. However, the dynamics of 
acoustic waves will be different from those of Schrodinger waves due to the frequency dependence of u. 
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B. The t matrix and resonant point scatterers 



To elucidate the notion of the t matrix, we start the microscopic description with a simple example: scattering of 
an incoming beam from one scatterer in one dimension. The quantum-mechanical wave equation in d = I with a 
scatterer in x = xq reads 

- ^"{x) - u6{x - xo)'i'ix) = E-ii{x) with E = . (74) 
Assuming that a wave e*'"'^ comes in from a; = — oo, we look for a solution of the form 

The constants A and B are determined by the continuity of ij) and the cusp in ij)' bX x — xq. The solution can then 
be represented as 

u + 2ik 

= e*'=" + G(2;,a;o)ie*'="°. (76) 

The first term is just the incoming wave. The second term represents the wave at the scatterer in xq, where it picks 
up a scattering factor t, 

l-m/(2fc)' ^^^^ 

and is transported to x by the medium without scatterers. The Green's function in the medium without scatterers, 
defined as {-d"^ /dx^ - P)G{x, x') = 5{x - x'), reads 

G(x,^')^^^. (78) 

Note that we have chosen the sign of ie such that the Fourier representation of the pure Green's function reads 
G{p) = l/ip'^P-ie). 



C. Tlie t matrix as a series of returns 

In the above example Eqs. ( [77| ) and (|7^) show that the t matrix can be expressed as 

i(^o) - z ^, r. (79) 

1 - uG{Xo,Xq) 

Expansion in powers of u yields the "Born series" , a series that have a clear physical interpretation: the term of order 
describes waves that arrive at the scatterer, return to it fc — 1 times, and then leave it for good: 

t= u + uGu +uGuGu + uGuGuGu + ... (80) 
• = X + X — X + >^^^ + X X X x + ... (81) 

The curved lines in the figure indicate that scattering occurs from the same scatterer (so far there is only one). The 
t matrix is indicated by a circle in the figure. For an extended scatterer this expression reads 

t(r,r') = -V{r)6{r - r') + V{r)G{r,r')V{r') 

d3r"y(r)G(r, r")V{r")G{r" ,r')V{r') + ... (82) 



Except for the case of point scatterers, this iterative solution of the single scatterer problem is not very helpful. The 
problem of vector wave scattering from a sphere was solved by Mie (1908), but the result is quite involved, too involved 
for our multiple scattering purposes. The approach shows, however, that the t matrix depends on the Green's function 
G, which describes propagation in the medium without that scatterer. However, the Green's function will depend 
on further details of the scattering medium, such as the presence of walls or other scatterers. Generally, it cannot 
be taken from the literature but has to be calculated for the problem under consideration: the t matrix describes 
scattering in a local environment. The standard infinite-space t matrix found in the literature only applies to that 
specific situation. 
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D. Point scatterer in three dimensions 

In three dimensions the Green's function reads 

G{r, r') ^ ' '/I - . I ^ n +ir + - r'|), (83) 

47r|r-r'| 47r|r-r'| in 

which is closely related to the Yukawa potential for hadron-hadron interactions. The divergence at r' = r will cause 
some problems in taking the point limit. 

1. Second-order Born approximation 

For weak scatterers one often truncates the Born series. A complication is that the Green's function in three 
dimensions diverges for r — > r', which makes the Born series (|8^) ill defined for point scatterers. In reality these 
divergences are cut off by the physical size of the scatterer, so they play no role for weak scattering. One thus keeps 
the first-order term and the imaginary part of the second-order term; the regularized real part of the second-order 
Born term will be small compared to the first-order term. This leads to the second-order Born approximation 

t = u + iuHmG{r,r). (84) 

For the system under consideration this becomes 

k 

t = u + iu^ — . (85) 



The fact that Imi > signifies that t matrix still describes scattering [see Eq. (106)]. More detailed aspects of the 
scatterer, such as resonances, are not taken into account. For electrons the second order Born approximation is often 
applied. Instead of working with point scatterers one usually considers a Gaussian random potential with average 
zero {V{y)) — and correlation {V{y)V{y')) = u^S{r — r'). This leads exactly to t = iu^ lmG{r,r)—iu^k/4:TT. For 
light, however, the second-order Born approximation is less applicable if resonances occur. 



2. Regularization of the return Green's function 



For finite-size scatterers the divergence of the Green's function at coinciding points is not a severe problem. The 
physical scatterer always has some finite radius, 1/A, that cuts off any divergency. When we wish to consider 
point scatterers, the divergence does cause a mathematical problem. In contrast to high energy physics, cutoffs in 
condensed matter physics represent physical parameters. Indeed, we shall identify A as an internal parameter of the 
point scatterer, that fixes the resonance frequency. 

The return Green's function G(r r') resulting from (p3) diverges. It should be regularized as 



G(r,r') 



1 



A- 



ik 
in' 



47r|r — r' 

Another method, often used in quantum field theory, is to introduce a large momentum cutoff, 

d^t 



(86) 



d^p 1 



G'''S{r,r) = A 



(27r)3 p2 

d3p 



+ 



^{G(p) - -2 } 



(2n) 



P 



1 



(o .,{G{p)--} = A 



ik 
in' 



(87) 



This subtraction is also possible if there are other scatterers or walls. Different regularization schemes define different 
point scatterers. The physical result is largely insensitive to such details. 
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3. Resonances 



If we insert the regularized return Green's funetion in Eq. (|79|), we ean write the t matrix in the same form as in 
the d=l case, by introducing the effective scattering length Ucs, 

II. 

t = 



1 - ■uG'™s(r,r) 
u 

l-"(A+f ) 

f/eft 



-IT: (88) 

1 7-7- ifc ' ^ 

where 

C/eff ^ (89) 
1 — uA 

For Schrodinger waves, the regularization brings just a shift in u as compared to the second-order Born approximation. 
This has no further consequences, since u and A are constants, and so is C/eff- For scalar waves there is an important 
difference, since then u is frequency dependent [c.f. Eq. (^3|)] : u = ak^. Let us call k^, the wave number at resonance. 
If we identifylj A = l/(a/c^), we get 



(90) 



One sees that UcS — *■ 00 for /c ^ /c*. Then t becomes equal to ~ Ani/k — 2iX. As A ^ 1/A this shows resonance 
with strong scattering: the effective scattering length \t\ = 2X is much larger than the physical size of the scatterer 
^1/A. 

The resonance is an internal resonance of the scatterer, comparable with the s-resonance of a Mie sphere. It is 
strongly influenced by the environment in which the sphere is embedded. For small frequencies t « q/c^, leading to 
the Rayleigh law ct ~ for w — > 0. 

4. Comparison with iViie scattering for scalar waves 

We compare the above result of regularization with an exact result. The t matrix for scalar s-wave scattering reads 
(see, e.g. Merzbacher (1970), page 238). 

47re~^*'^° 47re~*'^° sinfca , , 

f — (92) 

mkcot{mka) — ik k 

where a is the radius of the sphere and m the ratio between refractive indices of the sphere and the outside. The first 
resonance occurs at wave number fc^ = Tr/(2ma). It becomes sharp if one takes small a, large m, such that k^, remains 
fixed. For k close to fc* one gets 



-1 



For our point scatterer (^ we have 

1 - iUesk/4TT ^ ^" y a^k^ kl 



t= . -4.fl!^(l-l)-.fc) , (94) 



^That can be done provided a > 0, that is to say, when Escattcror > Emcdiun 
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where the scattering length is 



dV(e(r) - 1) = / d?r[m^ - 1] - ^(m^ - 1). (95) 



In the limit of small a this becomes a = '^^"g™ . Comparing Eq. ( |9^ ) with Eq. ( |9^ ) we find the prefactors At: /a and 
7r*/32m^a'^. The difference is a factor 7r'*/96 « 1.0147, thus the results coincide within 2%. In Eq. ( ^ we called 1/A 
a measure of the radius of the scatterer. We find, if m » 1, 

1 fc?(m2 - l)a3 TT^ 

- = -^^ — w — a w 0.822a. (96) 

A 3 12 ^ ' 

Indeed 1 /A is a good measure of the radius. We have thus found a simple expression for the t matrix which incorporates 
the essential physics near the s-wave resonance. Finally, we mention that Mie scatterers absorb radiation if the 
refractive index has a small imaginary part. For point scatterers absorption is described by giving u a small negative 
imaginary part. 



E. Cross sections and the albedo 



In the literature several cross sections are encountered. Here we discuss the three most important. Let a plane 
wave be incident in direction n'. As mentioned above, the total wave scattered from one scatterer is 



^'(r) = e^'^'"' '" + j (fY'dPr"G{r, Y')t{r', r")e*'="' '"", (97) 
with a real-space representation of the t matrix [see Eq. ( ^ . Its Fourier transform is called the off-shell t matrix, 

i(p,p') = J d^rd^r'e-'P"+'P' '''t{r,r'). (98) 



Let us assume that the center of the scatterer is located at the origin. Far away it holds that 

^ikr—ikn-v' 



G(r,r') 



r 



We insert this in Eq. ( |97| ) and with ( |9q ) we find that the scattered wave has the form 



^scij) « —t{kn,kn'). 



(99) 



(100) 



The scattering cross section is defined as the scattered intensity integrated over a sphere, normalized by the incoming 
intensity 



47r 



t{kn, fcn') 



Airr 



1 



(47r) 



dn\t{kn, fcn')|^ 



(101) 



4Tr 



Notice that the t matrix is only needed for momenta |p| = fc ("far field", "on the mass shell", "on shell t matrix"). 
One often denotes it by t{n,n'). An isotropic point scatterer thus has scattering cross section 



tt 
An' 



(102) 



The second important quantity is the extinction cross section. It tells how much intensity is lost from the incoming 
beam. Let us assume that a plane wave is incident along the z axis from z = —oo. We consider the intensity in a 
small solid angle around the z axis for large positive z. Because then r fa z + (x^ + y'^)/2z, the total wave is given by 



4-(r) = e'"' + tin, n)— e'^^^^ +y 



(103) 



The intensity for large z is 
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= ! + Re^^^e''=(^'+^')/2^ (104) 
2ttz 

Wc integrate this over x and y inside an area A perpendicular to the z-as. The condition that x and y be much 
smaller than z makes the integrals Gaussian, after which the z dependence disappears: 

/ dxdy'S*^ = A - aex- (105) 
J A 

The extinction cross section is the surface over which the incoming beam has to be integrated to collect an equal 
amount of intensity. It is equal to 

, . lmt(n,n) , 

fTex(n) = (106) 

For a point scatterer one has t{n, n) = t. The albedo of the scatterer is the ratio of scattered and extinct intensity, 

a=^ = -^. (107) 
aex 47rlmt 

For pure scattering a = 1; this is called the optical theorem. If absorption is present the absorption cross section can 
be defined as aabs = Cex — (^sc The albedo then equals 

a= ^ . (108) 

O'sc ~t~ ^ahs 

For an extended spherical scatterer the (extinction cross section is angle dependent. One defines the earlier encoun- 
tered phase function in terms of the t matrix as 

Kcose)=p(n.„') = ^^^, (109) 

with 

^p(n,n')=a. (110) 



/ 



Note that, because of the far-field construction, the optical theorem cannot be applied immediately in a system with 
many scatterers. Instead one has to impose the Ward identity, which is its generalization. 



V. GREEN'S FUNCTIONS IN DISORDERED SYSTEMS 

After our microscopic treatment of a single scatterer, we now consider scattering from many scatterers. The Green's 
function of a given sample depends on the realization of disorder: the location and the orientation of scatterers. Av- 
eraged over disorder it is called the amplitude Green's function. It describes, on the average, unscattered propagation, 
such as an incoming beam or the wave scattered from any given scatterer. It should be contrasted with the diffuse 
intensity, which is the (multiple) scattered intensity that will be discussed below. For a introduction to Green's 
functions in disordered systems see Economou (1990). 



A. Diagrammatic expansion of the self-energy 

The amplitude Green's function is related to two important concepts: the density of states and the self-energy. Let 
us c:onsider waves with frequency ^' = cko- It will be seen that the presence of many scatterers will change the "bare" 
wavenumber Uq into the "effective wavenumber" 

K = k + ^. (Ill) 

This describes a phase velocity Vph = uj/k = cko/k . 

For an electron in a random potential the definition of the amplitude Green's function is 
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Gr.r' =n ydK"P(^r")5r,r', 



with 



(112) 



(113) 



r.r 



The effect of the random potential is to introduce in the bare Green's function g{p) the self-energy S, 

^^^^ ^ P^-E-npy ^'''^ 

Exact averaging over the disorder is only possible in particular cases in one dimension and cannot be done in general. 
□ Therefore, we employ a diagrammatic approximation to calculate the self-energy. Suppose that we have a system 
with randomly located strong scatterers (water drops in fog, lipid particles in milk, Ti02 particles in a liquid or solid 
sample). In the limit of point scatterers the potential becomes V{y) = — uS{r — Hi). The Green's function of 
this problem is represented diagrammatically in Fig. ^. 

Go denotes a bare propagator, that is to say, the propagator in the medium without scatterers. The quantity 
needed is g, the Green's function of the random medium, also called dressed propagator. To calculate all diagrams 
would amount to solving the problem exactly, which is usually not possible. We shall therefore assume that the 
density of scatterers n is small. This will allow us to set up a perturbative expansion of S. We average over all 
disorder configurations, that is to say, all positions (and orientations) of the N scatterers in the volume V. We sort 
the diagrams (see Fig. ^) and get 

G = (g) = Go + GqYiGq + GoSGqSGo -)- . • . 

= Go + GoSG (115) 



1 1 



Go' - S -V2 - fc2 



(116) 



The relation ( 115 ) is the Dyson equation. The lowest order approximation to S was already calculated in Sec. IV, 
as it is proportional to the t matrix. Due to the averaging over the scatterer positions, there is a also a density 
dependence, which can already be calculated from the first order Born term. 



fd^K, d^RN ^ p. ^ fd^K uN 



(117) 



The same argument holds for all orders in the Born series, so one obtains to 0{n) 

Y.^nt. (118) 

This lowest order approximation is also known as the independent- scatterer approximation. In Fig. |^ it corresponds 
to the approximation E = •. The effective wave number is extracted by comparing the dressed and bare propagators 



K=^kl + ^^k+^^, (119) 



which leads to the complex index of refraction 




(120) 



^Exactly solvable models exist in one dimension. B.I. Halperin considers Gaussian white noise potential on a line (Halperin, 
1965); Th.M. Nieuwenhuizen considers exponential distributions on a ID lattice (Nieuwenhuizen, 1983; Nieuwenhuizen, 1984). 
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B. Self-consistency 

It is impossible to calculate all diagrams. Often it is feasible, however, to take into account, without much extra 
effort, all higher order contributions of a certain type. In the case of the self-energy we can take t = u/[l — uG(r, r)], 
instead of t = — uGo{r, r)]. This is called a self-consistent approach. Now the last diagram in the series for S of 
Fig. ^ is included. Similar terms with any number of intermediate dots are also accounted for. (One has to be careful 
to avoid overcounting of diagrams, however). Nevertheless, some two-scatterer diagrams have still been neglected at 
this self-consistent one-scatterer level. 

Physically the self-consistent method is very natural: it describes that propagation from one scatterer to another 
does not happen in empty space, but in a space filled with other scatterers. Therefore self-consistency is a fundamental 
concept, and it will serve to satisfy conservation laws (Ward identities). Indeed, without self-consistency there is no 
exact cancellation in the Ward identities. For instance, in the second-order Born approximation one will derive 
Eq. ( |128| ) with 1 — a ~ u*^, formally describing absorption or even creation of intensity, in situations where intensity 
should be strictly conserved. In order to work with such approaches one must neglect terms. In the more physical 
self-consistent approach such ad hoc manipulations are not needed and not allowed. 

The self-consistent t matrix in the independent scatterer approximation. 



1 - uA - j^^uy/k^ +nt 

can lead to a real value of t. This amounts in this approximation to a gap in the density of states (Polishchuk, Burin 
and Maksimov, 1990). Probably there is no real gap, but still a small density of states, which may lead to Anderson 
localization (Polishchuk et ai, 1990). 



VI. TRANSPORT IN INFINITE MEDIA: ISOTROPIC SCATTERING 

Just as the amplitude Green's function follows from solving the Dyson equation, the intensity follows from solving 
the Bethe-Salpeter equation. In this section we restrict ourselves to an approximate form of this equation: the so- 
called ladder approximation. It will be seen that this is corresponds to the independent scatterer approximation at the 
intensity level. The ladder approximation allows a microscopic derivation of the radiative transfer equation discussed 
in Sec. 



Ill 



In infinite media there are no boundaries which simplifies the analysis. In practice this situation applies to cases 
where the distance to boundaries is many mean free paths, so in the bulk of the multiple scattering medium. We first 
consider this situation and derive the diffusion equation from microscopies. In doing so, we find expressions for the 
diffusion coefficient and the speed of transport. 



A. Ladder approximation to the Bethe-Salpeter equation 

For the transport of energy we must consider the intensity. We have to multiply the Green's function by its complex 
conjugate. In other words, we must multiply the retarded Green's function by the advanced one. This must be done 
before averaging over disorder. We indicate this in Fig. |5| by drawing the expansion for g (first line) and drawing the 
one for g* below it. 

When scattering from a certain scatterer occurs more than once, we indicate this again by dashed lines. Visits 
to scatterers are included in g or g* , as before, but it may also happen that both g and g* scatter from a common 
scatterer. This leads to the connection lines between the upper and lower propagators in the figure. 

Of special importance are the ladder diagrams, depicted in Fig. |5|(a), which describe the diffuse intensity in the 
independent scatterer approximation. They lead to the classical picture of propagation of intensity from one common 
scatterer to another. As the intermediate propagation takes place in a medium with many other scatterers, each 
intermediate g or g* line visits many new scatterers: they are the dressed propagators of the previous section. 

Alternative names for the ladder diagrams are: ladder sum, diffuson, particle-hole channel. They are constructed 



in three steps: (i) Sum all immediate returns to the scatterers. As was explained in Sec. IV, this replaces the bare 
scattering potential u by the t matrix, (ii) Keep only those diagrams in which g and g* visit a certain series of common 
scatterers once and only once and in the same sequence. The intensity g*g propagates from one scatterer to another, 
(iii) Use the dressed Green's function as intermediate propagator. As a result in between the common scatterers both 
g and g* visit any amount of other scatterers. 
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Let /(r) = |4'(r)p denote the intensity arriving at point r. According to Fig. || it can be decomposed into terms 
without scattering, terms with one common scattering, term with two, etc: 

/(r) = |vl/„,(r)|2+ntt J dV|G(r - r')n*™(r')P 

+ intt)2 J d'v'd^v"\G{v ~ v')\^\G{v' - r")n*™(r")P + • • • (122) 

which can be written as an integral equation 

C{v) = ntt\^,n{Y)\^ +ntt j (fr'G{r-v')G*{r-v')C{r'). (123) 

where £(r) = nttl{r). The ladder propagator £(r) is the intensity that leaves point r after being scattered. 

In the second-order Born approximation the same scatterer is at most visited twice consecutively, with the result 
that the t matrices in the above ladder equations are replaced by the scattering potentials u. 

B. Diffusion from the stationary ladder equation 



In the bulk, that is to say far from boundaries, the source term in the ladder equation vanishes. Eq. (123) therefore 
takes the form 



£(r) = ntt I I G(r') 1^ C{r + v')dPv' . (124) 



Now assume that C{v) varies slowly on the scale of one mean free path and expand 



£(r' + r) = C{y) + r' • \I C{r) + ^rV : VV£(r). (125) 



where : denotes a tensor contraction. Inserting this in Eq. (124) yields three contributions. The first term is 



-r'/l _ £ 



nttC{v) / (Py'- — -- ^ ntt—C{Y) = aClr). (126) 
J {'iirr'Y 47r 

In the prefactor we have recognized the albedo a (see Sec. IV). The second term vanishes due to symmetry. The third 
term yields 

inft j (fr'G*{r')G{r')Y'v' : VV/:(r) = ^ntt x i x 2e'^'^C{r). (127) 



Inserting this in Eq. (124) yields the stationary diffusion equation 



V^£(r) = ^-^^C{v). (128) 
in agreement with eqs. ( p^ and (|57|), since for isotropic scattering it holds that l^c — ^tr- 

C. Diffusion coefficient and the speed of transport 

Here we give a simple derivation of the nonstationary diffusion equation. In doing so, we automatically encounter 
the speed of transport. We again consider isotropic scatterers. In an infinite system it is useful to consider the 
Fourier-Laplace transformed ladder equation 



£(q,f]) = 5(q,r!)+nt>_)t(^+) j ^G(p+,c.+ )G*(p_,c._)£(q,r!) 

^(q, 17) + A./(q, 0)/:(q, O) 
5(q,r!) 



1-M(q,.)' (^2^) 
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where S is the transformed source term; its precise form is of no interest here. We introduced 



(130) 



The parameters p and uj are "internal" or "fast" variables, which involve one period of the wave, q and fl are 
macroscopic or slowly varying parameters. They describe variations over distances much larger than the wavelength 
and times much larger than the oscillation period. 
At = 0, the bulk kernel has the property 



M(q, 17 = 0) = ntt 



d^P ^/ , \r'*f \ arctan(g£) 
G(p + q)G (p) = 



(2^) 



qi 



(131) 



We want to know the behavior for large distances and times, so for small q and il. We expand to orders H. and q^. 
Denoting G = G'(p, lj), it holds that 



G 



9 Iq^^^^^ dt\ 



+ (p • q)2G3 



Inserting this in Eq. (129) implies for the kernel 



, Ll> nt' 



uj nt' , 



/n - n{- + — )/2i + n{- + —)h 



^2 ^2^2 
--^(Il2 + -^21) H ^ — (-^31 + A3 ^ I22) 



Here we have defined the integrals 



kl 



d'p G^-(p)G*'(p) 



(27r) 



They are calculated in the appendix. Inserting their values yields 

M(q, n)=a- iilTsc - -q^f. 

3 

As discussed by van Albada et al. (1991), r has two contributions, 

~ H~ '^dw — '^sc '^dw • 



(132) 



(133) 



(134) 



(135) 



(136) 



The first term is the "scattering time" or "time of travel" — ^1 c. The explicit expression for the "dwell time" t^iu 
follows as 



^ dlogi 27r ^ dt 
Tdw = Im— h -r^Re— . 

duj kott duj 

We insert the i-matrix of a point scatterer, Eq. (^), yielding 



(137) 



incurs 47r[/2 



off 



ko 



kott Vf^eff 



kott 



47r 



koUcff 



duj 



(138) 



Notice that the terms without i/^g have compensated each other. This cancellation follows more generally from a 
Ward identity (van Albada et al., 1991). For electrons in a random potential Tdw therefore vanishes. 
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For acoustic waves and light waves the situation is more interesting. Since UeS depends expHcitly on frequency, C/^g 
does not vanish. Both terms are additive, leading to a finite dwell time Tdw- Using Eq. this yields at resonance 

2tt 1 period 

Tdw = — • —TT = r. — ■ (139 

u) ak"^ coupling 

This dwell time becomes longer when the coupling to the environment, that is, the normalized scattering strength 
k^a — 47rfc'^a'^(m^ — 1)/3, becomes weaker. Therefore this is an important experimental effect. The speed of transport, 

(140) 



can be substantially smaller than the speed of light when realistic values of n are inserted in this formula. 

For resonant atoms at fixed positions, the reduction of speed may be as large as 10^ (Vdovin and Galitskii, 
1967; Nieuwenhuizen, Burin, Kagan and Shlyapnikov, 1994). As a final result the diffusion coefficient D = vi/3 
reduces, as is observed experimentally. For the general situation and for more details, see van Albada and Lagendijk 
(1985) and Lagendijk and van Tiggelen (1996). 

With this result for t we find for the ladder propagator ( |129| ) at small q, fl 

r«,,o).aa2M^_^^. ,141, 



involving the reduced external frequency 



n^^. (142) 



This is exactly the propagator of Sec. || with absorption length Labs = 1/k = — a). 

One often considers the Schwarzschild-Milne equation with stationary S source, S{r,t)— nttS{r)^4:Tr£~^S{r). In the 
diffusion approximation this yields 

Cici,n)^^- (143) 

This form is commonly called the "diffuson." 



VII. TRANSPORT IN A SEMI-INFINITE MEDIUM 

In this section we consider in detail the transport equation for a very thick slab. We follow the approach of 
Nieuwenhuizen and Luck (1993). We shall also consider the case of a mismatch in refractive index between the 
medium and the outside. 



A. Plane wave incident on a semi-infinite medium 



In many situations the index of refraction of the scattering medium differs from that of its surroundings. An 
example is Ti02 particles suspended in a liquid. Boundaries between media of different indexes act partly as mirrors. 
They cause a direct reflection of the incoming light and re-inject part of the multiple scattered light that tries to exit 
the system. As long as the system is optically thick, they lead to effects of order unity, no matter how small the ratio 
X/£, thus they are important in a quantitative analysis. 

We shall calculate the angular resolved intensity profile for the case of a plane wave incident on a perfectly flat 
interface. This will create a specular reflection at the interface. The situation of nonspecular reflections from nonideal 
interfaces also has some practical relevance. 

Consider a semi- infinite medium. For z > there is a scattering medium with refractive index no. For z < there 
is a dielectric with index m = no/m. Our notation is indicated in Fig. |l|. The system is governed by Eq. (73). We 
first determine the incoming wave in the scattering medium. We consider a system with scatterers in the half-space 
z > 0. We replace the action of the scatterers by a self-energy term. The average wave equation reads, after Fourier 
transformation of the transverse vector to q_L = {qx,<ly), 
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— ■^{z) + PH{z) = Q, z>0 (144) 
(P 

— ^'(z)+p2^-(z) = 0, z<0 (145) 
P^^k^-q±^ + nt, p^^kl-qj_^. (146) 

A plane wave of unit amplitude, incident from z < 0, causes a reflected wave for z < and a refracted wave for z > 0: 

Pa+Pa 

g.kl.p+»p„^ (z>0). (147) 



Pa+Pa 



The prefactors in Eq. (147) follow from the requirement of continuity of ^E* and its derivative. To lowest order in the 
density we find for the real and imaginary parts of P 

P^fcocos.' + .^, (148) 

where 9' is the angle of the refracted incoming wave with respect to the z axis. 

The source term in the scattering medium, the unscattered incoming intensity, can be written as 

/,„ = |vi/,„|2 =1 |2 e-/^cos.„ ^ E2.T,e-^/i^^. (149) 

-^a ' Pa ^a 

Inside the random medium the unscattered intensity is damped exponentially (Lambert-Beer law). To leading order 
in 1/kl, 



11-^11^-1 + 1/r 



(150) 



yu + ^p? -l + l/n, 

where R{^) is the angular reflection coefficient for scalar waves, and R — \ for total reflection, which occurs when 
m > 1/ yl — ^j? . For vector waves the same equation applies in the s-wave channel (Amic et al., 1997). In this 
expression for R{fJ.) the imaginary part of Pa, of order l/k£, has been neglected. It has been shown by Nieuwenhuizen 
and Luck (1993) that we must neglect this to ensure flux conservation. 

We use the optical depth t = zjl (not to be confused with the average time per scattering r = + Tdw) and 
introduce r(r) as 

I[z)^^L{z^Tl) = ^TanT). (151) 

F is the diffuse intensity per unit of intensity entering inside the medium. In terms of F the ladder equation becomes 
dimensionless 

/■OO 

F(r) = e-^/^° + / dT'A/(T,r')r(r'), (152) 



where M(r, r') follows from the square of the amplitude Green's function. In contrast to Eq. ( |36| ) G now consists of 
two terms for z, z' > 0: a direct term and a term involving reflection for the boundary z = 0: 

G(z, z', q^) = ^{e^^l^-^'l + ^e'^(^+-')}. (153) 

We shall later need G for z < 0, z' > for which 

G(z,z',qi)--^e-'f^+^^^'. (154) 
P + V 

Inside the random medium, | G P thus has four terms. The cross terms oscillate quickly and can be omitted. We 
thus keep the previous bulk contribution to the kernel Mb and the new layer kernel M^, M — Mb + Mj^. For the 
bulk kernel we find 
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MBiz,z') = 47r / d^p I G(z,p;z',p') p . 



It holds that 



Since 



one gets 



Mb{z,0) = / dxdy 



2 ^-^x^+y^+zye 



z + p = z /fi => pdp = -z d/i//r 



(155) 



(156) 



(157) 



Ms(z,0) = 4^ [dxdv—^e~'-/'^ = / 



2/i 



(47rr)2 



pdp 
2r2 



yieldmg for the bulk kernel [cf, Eq.(p6[)] 

Mb{t, r') = £ |^c-l--'l/^ = i£;i(|r - r'|), 



3 2/i 

where E\ is an exponential integral. In the same fashion the layer term is 

Ml{t,t') = 47r / dxdy 



(158) 



(159) 



i P — p 



P{z+z') 



2p 



2PP + P 



(160) 



It consists of three effects: exponential decay of intensity that goes from the point r' with angle 9 towards the wall 
having an optical path length t' / p,; a reflection factor i?(/i); a fm'thcr decay over an optical path length r/p between 
the wall and the observation point r. 

We now calculate the angular resolved reflection. For a plane wave incidence the solution of Eq. ( |123| ) is C{z, p) — 
C{z). We can write for the reflected intensity at point (z', p') with z' > 



Ir{z\p') = / dzd^p£(z,p) I G{z,p;z',p ) f 



(161) 



Here C{z) is related to T{t) via z = £t and Eq. ( [1511) . 

It is useful to consider first a finite area A and to continue this area periodically. We can then integrate out the p 
dependence; later we shall take the limit A oo. Using the periodicity of the surface one can express the Green's 
function as 



■ip-p') 



(162) 



Since the incoming plane wave is infinitely broad, the diffuse intensity I does not depend on p. The integration over 
p is simple: 



I d'p I G(z, p; z', p') 1 ^ GG* « I G(z, z', q^) | 



(163) 



The last step holds in the limit A^ oo. If we denote the outgoing qx by q*]^, we can express d^qj_ as 

d^qi = q\_<^h\_<^^b = kj sin Ob cos 6lf,d6lbd0t = kj cos ebdfli,. (164) 



Inserting Eqs. (164) and (|l63|) in Eq. (|l6l|) yields a result that is independent of the observation point (z',p') 
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In(Ap') = /dz/ ^^df^^^:r^/:(-)e-^/^™^^ (165) 

where 0^ is the dhection of the radiation in the medium, that is refracted into the outgoing direction 0b [see Fig. 
we also recall that fib = cos 6*1,. We finally find for the angular resolved reflected diffuse intensity 

^ = ^ ^ ^^JL [ dz£(z)e-/^™^^^ (166) 

dilb ailb 47r-^ PbPb J 



Using Eqs. (151) and ( |146| ) at nt — > 0, one can express the numerical prefactor of the integral 

■ cos 0b4:Tr—- — = cos r 



APbPb (2^)2 AT:pb "PaPb 

COS0,— -— = — ^— -. (167) 



47rcos6'6 ko^ako^ib Airm'^ fiafJ-b 

For the angular resolved diffuse reflection of a semi-infinite medium we thus find 

cos TO 
47rm2 yUaMb 

with the generalized bistatic coefficient 



ABidaM = TZZ^-^lif^a, Aib), (168) 



7(Ma,Mb)=/ dTrs{fia,T)e-^^^' = drdr' G s{r , T')e-^ ' b'^' ' . (169) 
Jo Jo 



In the absence of index mismatch (jn = 1) one has Ta^b = I7 ^■a,b = cos^a^b, so the prefactor in Eq. ( |168| ) becomes 
l/(47rcos0;,). Fig. ^ shows numerical results for Aji for perpendicular incidence {9a = 0), see (Nieuwenhuizen and 
Luck, 1993). 



B. Air-glass-medium interface 

We now consider a semi-infinite scattering medium, separated from the air by another dielectric, such as glass, of 
thickness d. For z < —d there is air, for —d < z < glass, and for z > the scattering medium. A plane wave of 
unit amplitude comes in from z = —00 with wave vector (qj_,fcz). The perpendicular component q_L is conserved 
at the interfaces. We denote by pi the z-component of the wave vectors in the three sectors i = 0, 1, 2, where i = 
corresponds to the scattering medium, i = 1 to air, and i = 2 to glass. The wave numbers in the three media are 
ko ~ oj/cq, kl — Lo I c\ ,and ^2 — oj / C2, respectively, where Ci is the speed of propagation in the medium i. The 
incoming, refracted, and specularly reflected waves are given by 

^'(r) = { tie'^i^-P+'P^^ + ne'^i^-P-'P^^ (-d < z < 0) 
te^q^-p+voz (2 > 0). 



Here r is the reflection amplitude of the system, t the transmission amplitude, and pi — \/A;| — qj.^. Continuity 
requirements at z — — d and z — yield r, ri and t, ti: 



r 



{Po + P2){pi - P2) - {po-P2)ipi +P2)e''P^'' 

{P0+P2){Pl +P2) - {P0-P2){Pl -p2)e2^P2d 



^ ^ 4piP2e'(P^-P^)'^ 

(PO + P2)(P1 + P2) - (Po - P2) (Pi - P2)e2'P^d ■ 



The reflection and transmission coefficients of the double interface are 

|rp ; ^|tp, (171) 
Pi 

respectively. As we assume the thickness of the glass plate is not smooth within one wavelength, these expressions 
must be averaged over the spread in thickness. This amounts to averaging over the phase 2p2d = ip and leads to the 
average transmission and reflection coefficients 
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Po f dip ^ ^ 2 



, ,-, • (172) 

The integrand is of the form A/{B — Ccos</)). The final result is 

n-i-R.- ^ wip. (173) 

[Pa + Pi){PqPi +Pi) 

We can check this in special cases. Inserting P2 = Pi — P and po = P indeed reduces to previous the result for two 
media. 

We can describe this system by replacing T(/i) in previous equations by T^{pi) and replacing the reflection coefficient 
_R(/i) in the Milne- kernel by R-i{^). As before, ^ = po/ko is the cosine of the angle 6' between the radiation and the 
z axis. 

Specular reflections in the glass have now been taken into account. For a not too narrow beam this is useful for 
thin glass plates. For a not so broad beam impinging on a medium with thick plates, multiple reflections from the 
glass interfaces can results in components that fall outside the incoming beam and even outside the medium (Ospeck 
and Fraden, 1994). 



C. Solutions of the Schwarzschild-Milne equation 

We consider properties of the transport equation in a semi- infinite space (Nieuwenhuizen and Luck, 1993). The 
Milne equation (152) has a special solution Ts{pa',T), while the associated homogeneous equation without source 
term has a solution Th{t). We are interested in the asymptotic behavior of Ts{fJ.a'iT) and r//(r). Deep in the bulk 
{t ^ 1) one expects a slow variation of T{t') at the scale of one mean free path (|r — r'| =0(1)) and one can expand 

r(r') = r(r) + (r' - r)r'(r) + i(r' - r)2r"(r) + • • • (174) 



When this is inserted in Eq. (152), one finds for r ^ 1 

4- 
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r(r)«r(T) + ir"(r) + .-- (175) 



Here we used that for large r 



The term 0{T') vanishes due to symmetry of the (r' — r)-integral. Eq. ( |l75| ) thus yields again the diffusion behavior 
T"{t) ~ 0. We therefore consider the homogeneous solution F^f and the special solution Ts with the asymptotic 
behaviors 

^h{t) ^ T + Tq . , , , 

J- ( \ ( \ V ^ °°)- (177) 

rs(T;Aia) ~ Ti(Ma) 

Corrections occur due to the interface and decay exponentially in r. For a numerical solution one may introduce 
(5r5(T) = r5(r) — T\{[ia) and JF/f (r) = V^ir) ~ t — tq and integrate them from r = oo. The requirement that they 
vanish there, determines tq and ri. The constant tq depends only on the index ratio to, while Ti(fj,a) also depends on 
the direction of the incoming beam. 

We define the Green's function Gs{t,t') as the solution of the inhomogencous equation 

/"OC 

GsiT,T')=6{T-T')+ dT"{MB{T,T")+ML{T,T")}Gs{T",T'), (178) 



subject to the boundary condition 6*5(00, t') < 00. This fmiction is symmetric, Gs{t,t') = Gs{t',t), and has the 
limit 

lim Gs(t,t') = ^Fh(t). (179) 
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The latter equality can be proven by taking the limit r' — > oo in Eq. (178). The delta function vanishes and the 
remaining equation is the same as the one for TniT). Therefore Gs{t,t' oo) is proportional to Th{t). The 
multiplicative prefactor can be fixed by expanding Gs{t,t') in t' — t and inserting this in the right-hand side of 
Eq. (178). Using Eq. (176) again, one finds for r, r' 1 



The solution is Gs{t,t') = 3min(T, r') in the regime (r, r', |t — r'| ^ 1). 
reads 1/3 in reduced units, that is to say, D = v£/3 in physical units. 
From Eq. ([l7^) it follows that 



(180) 



rs(T;/i) = / dT'Gs{T,T')e' 
Jo 



r' /n 



The diffusion coefficient D in Eq. (17£) 



(181) 



and, in particular, using eqs. (177) and (179), 

Ti{n) = lim rs(T;^) 



1 

D 



rH(T)c-^/^d7 



(182) 



The physical interpretation of Ti(/i) is the limit intensity {z — ri oo) of a semi-infinite medium. 

Numerical values of the injection depth tq, the normalized limit intensity ri (1) and the normalized bistatic coefficient 
7(1, 1) can be found in Table | for various values of the index ratio m. 

In Nieuwenhuizen and Luck (1993) it was verified explicitly that in this approach flux is conserved. This derivation 
will not be reproduced here. We refer the interested reader to the original paper. 



VIII. TRANSPORT THROUGH A SLAB 

We now discuss the transmission properties of an optically thick slab with isotropic scatterers. We derive the 
"ohmic" or diffusive scaling behavior T ^ £/L mentioned in Sec. ||, and give the full angular dependence of the 
transmission and reflection. This result is then used to calculate the resistance of an idealized conductor. 



A. Diffuse transmission 



We consider a medium with finite thickness L £. The medium has optical thickness b — L/£ ^ 1. For the 
moment we wish to neglect boundary effects. Therefore we are restricted to positions not too close to the boundary 
(10 mean free paths is a good measure, as the corrections decay exponentially). The solution for r(r) is a linear 
combination of the special and the homogeneous solutions, 



T{t) = rs(r) - aTniT) for < r < 10, 
= Ti(/z) — a{T + To) for r > 10. 

Near the other boundary it holds similarly that 

T{t) = a'Tnih - t) for < 5 - r < 10, 
= a{b - r + To) for 6 - T > 10. 

As they have the same functional form, both shapes can be matched in the bulk of the sample. This yields 

n(M) 



a = a = 



5 + 2to ' 



(183) 



(184) 



(185) 



Inserting this value in Eqs. ( 183 ) and (184) gives the intensity anywhere in the slab, expressed in terms of Ts 
and Th of the semi-infinite problem. Notice the important role played by the diffusive behavior in the bulk. The 
Schwarzschild-Milne equation has brought the precise behavior at scales of one mean free path from the boundaries. 
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To calculate the angular transmission profile, we follow the derivation for the diffuse reflection [see Eq. ( |168| . The 
expression for the differential transmission coefhcient per unit solid angle dQi, of a beam incident under angle 6a is 
given analogously as 

dr(«^ = SSl^f^ f drr(r)e-(^-)/- . (186) 
Using the solution for r(r) we rewrite the integral 



Jo J -oo 

/>oo 

— a dTTH{T)c 
Jo 



(b-T)//ib 



^JMD. (187) 



where we used Eq. (182). We have derived the angle-dependent differential transmission coefficient for a slab of optical 
thickness b — L/£, 

dTja^b) _ A^{eaJb) _ cosBaTgTb , , , . 

dnb ^ b + 2ro ^ 12nm^f,af^b{b + 2rof'^''''>^'^''''>- ^^^^^ 

As the intensity at the side of incidence (z = 0), equals r(r) = Tsir) — ar//(T), it is clear that the transmission term 
aTu{b — t) arises at the cost of the reflection. Therefore also flux conservation is also satisfied for an optically thick 
slab. 

In Fig. ^ we show A^{9a,0b)-, the normalized angular resolved transmission for perpendicular incidence {9a — 0), 
for several values of the index ratio m. 



B. Electrical conductance and contact resistance 

For metallic conductors in the mesoscopic regime the conductance is given by the Landauer formula (Landauer, 
1975) 

a.b 

where ft,/e^ ~ 25ki} is the quantum unit of resistance. This equation simply counts the weight of the channels that 
contribute to transmission. In the above description the wave number fcg can be replaced by the Fermi wave number 
vector kp. The analog of the contrast in refractive index is now the potential difference between the conductor and 
the contact regions (Vi 7^ Vq). In our formalism the conductance is given by 

2e^^ _2e2^ cos^_2e^_fc|^4_ 

^- h ^ ~ h 2^^'''cos9a~ h 37r(6 + 2ro)' ^ ^ 

with Tab = {Oa-, &b)/{b + 2to). Herewith one gets 

L + 2zo ' Stt/i ' 

where zq = tq£ and ctb is the "Boltzmann" value for the bulk conductivity. We distinguish a bulk resistance and a 
contact resistance i?c, 

The number of modes can be estimated as N Ri Akp. The above expression shows that Rc is proportional to the 
dimensionless thickness tq and inversely proportional to the number of channels. The nontrivial part is coded in tq. 
It depends on the potential drop V — Vi, but not on the density of scatterers. Because 
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kl-Vi = kl ,kl-Vo^kl (193) 



we can make an analogy with light scattering by putting kg = m kf 

, kl-v„ 



IX. THE ENHANCED BACKSCATTER CONE 

So far we have seen mainly effects that are largely diffusive and that could to some extent also be derived from 
particle diffusion. The enhanced backscatter cone is the clearest manifestation of interference due to the wave nature 
of the light. The wave character manifests itself most clearly in loop processes. The advanced and retarded waves 
can go around in two ways, in the same direction and in the opposite. This leads to an enhanced return to the origin, 
which is the basic mechanism for Anderson localization (Sec. ^. 

The optical enhanced backscatter in the exact backscatter direction has the same characteristics as a closed loop. 
It brings two possibilities, thus a factor of 2, for all scattering series involved. For Faraday active media it will be 
suppressed in a magnetic field. Away from the backscatter direction there is partial extinction. 



A. Milne kernel at nonzero transverse momentum 



We follow the discussion of Nieuwenhuizen and Luck (1993). The backscatter diagrams are closely related to the 
standard ladder diagrams. The only difference is the crossed attachment of the incoming and outgoing lines to the 
first and last scatterers. These diagrams are called maximally crossed diagrams (Fig. ||). 

Let the incoming wave be denoted by a and the outgoing wave by b. Their wave vectors have components 
(k^,pa), 0^±,Pb)- The product of incoming and outgoing waves at the first scatterer is 

Pa 

with p = {x,y). We define the transverse wave vector 

Q EE (kl - k^^). (196) 

For perpendicular incidence one has {9a = 0; k" =0). We consider the regime of angles close to the backscatter 
direction {9b — 0). Then it holds that |Q| ~ ki9b so that 

TO 

Consider the diffuse intensity / in the backscatter cone. It is given by 

/(r) ^ ffie''='-''e-^/^ + ^ /dr'|G(r - r')|'/(r'). (198) 

TO i J 

Notice that -dependence occurs only in the source term of this integral equation. Inserting /(r) = e^^'''I{z, Q) yields 
the Milne equation for I{z, Q), 

I{z, Q) = ^e-^/^ + / ^Mc(z, z', Q)I{z\ Q), (199) 

TO J £ 

with the Q-dependent Milne cone kernel 

Mc{z,z',Q)=A7T J d^p'e'^-'^p'-P^G{r,r')\^ . (200) 
The normalized intensity of the maximally crossed diagrams satisfies 
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mI{z,Q) 

rc(T, Q) = + y dr'Mc(T, r', Q)rc(T', Q). (201) 

The Milne kernel again has a bulk and a layer term, 

Mc{t,t',Q)^ MBiT,T',Q) + MLiT,T',Q). (202) 
Inserting |G'(r)p — e^'"/^/(47rr)-^ yields 

Mb(^,z',Q) = 4^ / pdpd0 .^ ^ »Qpcos0^-V(^-.o-+pV^ . (203) 

In terms of the variables t = z/i and /i = cos 9' , where ^' is the angle between p and the z axis, we find that 

^{z-z'y+p^ _ \t-t'\ pdp _ |T-r'|2 



-d^ . 



Inserting this in Eq. (203) yields 



Mb(|t-t'|,Q) 



Jo 2^7_, 2^ 

where Jo is the zeroth-order Bessel function. A similar analysis yields for the layer term 



(204) 



M^{t + r', Q) = 1^ |;^(m)Jo (Q^(t + t')./JF^^ c-^r+r')/,. (205) 



This form is useful at small and large Q. 



B. Shape of the backscatter cone 



The intensity in the backscatter direction consists of two parts: a diffuse background A^ , discussed in Sec. VII and 



a contribution A'^ from the maximally crossed diagrams. In the case of perpendicular incidence we have found for 



the background contribution Eq. (168) 



where the normalized bistatic coefficient 7 is given by 



nOO 

7(Ma,/^b)= / dTe--'^^T^{T,^la) (207) 
Jo 

) />oo 

dr / dr' e""/^" - Gg (r, r') . (208) 



/o Jo 

The contribution of the maximally crossed diagrams, normalized by the diffuse background, is given by 

A^iQ) = , (209) 

where 7tr is the amplitude of the paths that are identical to their time reversed analogs. Such paths do not yield 
a time-reversed contribution, and this should thus be subtracted from the first ter m. Fo r low scatterer density, 7tr 
consists of the single scattering event (in, scatter, out). This process yields [cf. Eq. ( |l69| )] 
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cso 



7TR = / dr / dr' c—-'5{t -t') = \. (210) 
The nontrivial term is defined as 

/•oo /-oo /-ex: 

7c(Q)-/ dTe-^Tc{T,Q)= dr / dr' e-^-^'Gc(r, r', Q) , (211) 
Jo Jo Jo 

where Gc is the solution of Eq.( ^Ol[ ) with source 6{t — r') instead of e^'^ and vanishing for t ^ oo. In Fig. |l^ we 
present the numerical results for A^' (Q) for several values of the index ratio m. 



C. Decay at large angles 



From Eqs. (204) and ( p05| ) it follows that Mb and Ml decay quickly as functions of r when Q is large. Physically 
this happens because t he e nhanced backscatter is suppressed at large angles due to dephasing. For large angles we can 
therefore restrict Eq. (211) to low order scattering. To second-order we have Gc{t, r'; Q) = 5{t — t') + Mc(t, t'; Q). 
This yields 



1 r°° 

7c(C?) = 2+ I drd 



r'e-^-^' 



For the bulk term 



f 

Jo 



dr dr' e 



Mb{\t-t'\;Q) 



{Mb{\t-t'\;Q) + Mi^{t + t';Q)) 
(Q) dfi d(j) 



(212) 



2^ 1 + l/;x - iQe cos (jjy/l/n'^ - 1 

^ d^ 1 



/o 2 ^(/i + l)2 + Q2^2(i_^2) 

Qlarge 7T 



(213) 



Inserting Eq. ( p05| ) yields for the layer term in Eq. ( |212| ) 

/•oo 

/ dr dT'e-^-^'ML(r + t';Q) 

JO 

'■Mm r 



:i?(A*)- 



TT 27r (Ai + 1 - iQi cos <^v/l-M^)^ 
A* + 1 



d/i fJ.R{p) 



[(Ai+l)2+g2^2(l_^2)]3/2- 



(214) 



The integrand decays quickly for large values of Q £ {1 — fi"). Setting fi — I ~ x/ (Q £ ) we obtain for the layer term 



Ril) 



dx 



This yields finally 



(4 + 2a;)3/2 2QH^ 

Tjl? ( ^ R{1) 
47rm2 \AQ£ 2Q^i^ 



(215) 



(216) 



for large Q. The effect of the layer term is essentially the square of the bulk term. This is because the paths involved 
are essentially twice as long (Fig. |ll|) 
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D. Behavior at small angles 

Assuming a linear Q dependence, we expand Tq around Q = 

rc(r; Q) ~ rc(T; Q = 0)- gf c(t; Q). 
Inserting this in Eq. ( |20l| ) and subtracting the (5 = terms yields, to leading order in Q, 



rc(T;Q) 



dr'A/c(T,T';Q)fc(T';Q) 



As this is a homogeneous equation, it holds that Vc{t\ Q) — aVYi{T; Q),so that Eq. ( pi7| ) becomes 

rc(r; Q) = rs(r) - aQrH(r) "i' ri(l) - aQ£r. 
Deep inside the scattering medium (r ^ 1) the diffusion approximation holds, implying 

r'^ir; Q) = QH^Tcir; Q) ^ Tc{t; Q) = ae"'?^^ . 



(217) 
(218) 

(219) 
(220) 



Matching this for small Qir with Eq. ( pi9[ ) one finds a = a = ''"i(l)- Using Eq. (207) this yields finally 

7c(Q)=7c(0)-3Q£r2(l) 



7c(0) 1 



Q_\ 
AQj 



where we used Eq. (182) with D = ^. The normalized opening angle AQ is 

^^-'^rf(f)- 



(221) 



(222) 



Note that 7c(0) = 7(Q = 0; /^a = Mb = 1)- The linearity of the peak of the cone is a result of diffusion, that is to say, 
of long light paths. The complicated expression for the opening angle shows that the skin layer plays an important 
quantitative role. The reason hereto is obvious: the light has to traverse it. 



X. EXACT SOLUTION OF THE SCHWARZSCHILD-MILNE EQUATION 

Some useful formulae were given in previous sections for calculation of the diffuse intensity. However, the 
Schwarzschild-Milne equation can be solved exactly in the absence of internal reflections (Chandrasekhar, 1960). 
We summarize the approach of Nieuwenhuizen and Luck (1993) here and generalize the Schwarzschild-Milne equation 
to include internal reflections. 



A. The homogeneous Milne equation 

The starting point of this analysis is the integral form of the radiative transfer equation, i.e., the Schwarzschild- 
Milne equation, Eq. ( 152 ). In the absence of internal reflections, = 0, the remaining kernel Mb depends only on 
the difference (r — r'), giving the Milne equation the structure of a convolution equation. Because of its half-space 
geometry the problem is still nontrivial. 

We con sider first the homogeneous Milne equation. Its solution Th{t) has the asymptotic behavior Th{t) = r + to 
[Eq. ( |177[ )]. We define the Laplace transforms of Mb(t, 0) and of Th{t) as 

/oo 
Mb{t,0) e'^'dr (-l<Re.s<l) (223) 
-OC 

9h{s)^ THir) e'^dT^DTi{^ = -l/s) (Res< 0). (224) 
Jo 

In the case of isotropic scattering the Milne kernel Eq. (|158|) leads to 
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The small- s behavior of m{s) 



1 1 + s 

m(s) = — In . 

^ ' 2s 1-s 



m{s) = 1 + Ds^ + 0{s'^) (s->0), 



(225) 



(226) 



determines the dimensionless diffusion constant D — 1/3. 

It turns out that the problem can be solved for any symmetric Milne kernel. The homogeneous Milne integral 
equation is equivalent to 



di_ m{t)gH{t) 

2'!Ti t — s 



(-1 < Ret < Res < 0), 



with 0(s) = 1 — to(s), and the asymptotic behavior Eq. ( |177| ) is equivalent to 

s^gH{s)-=l-Tos + 0{s^) (s^O). 



(227) 



(228) 



Equation (227) can be solved in closed form, for an arbitrary bulk kernel. Let us consider first the following "rational 
case," where Mb{t,0) is a finite superposition of iV decaying exponentials, namely. 



N 



AfB(r,0)=^ 



(229) 



with weights Wa > and decay rates (inverse decay lengths) pa > 0. We then have 



N 



a=l 



N 



(230) 



By comparing with Eq. ( P26[) , we obtain the normalization conditions ~ 1' ^a'^'^/Pa — Equation ( 227 ) 

can be evaluated by means of residue calculus. It loses its integral nature and becomes 



[s)gH{s) 



N 

E 

a=l 



2{pa + S) 



In order to solve Eq. (231), we first write the rational function (j}{s) in factorized form. 



^(s) 



Yltii^ -Pits' 



(231) 



(232) 



The {N — 1) zeros of |"(/i(s )/s^] in the variable s^ have been denoted by z^, with RcZq. > 0, for 1 < a < — 1. The 
normalization of Eq. (232) has been fixed by using <j){oo) = 1. 

An alternative expression for the diffusion constant D is derived by setting s zero in Eq. ( 232 ) and comparing with 
Eqs. dH) and (1230). We get 



D 



The solution of Eq. (|23l| ) is of the form 



4>{s)gH{s) = — — . 

Yla^iyPa + s) 



(233) 



(234) 



with A/'o(s) a polynomial of degree {N — 1), which can be determined as follows. Consider the value s — — Zq,: we 
have (/>(— Zq,) = 0, whereas gH{—Za) is finite. Equation (234) shows therefore that s = —Za is a zero of the polynomial 
Ao(s). Since there are exactly {N — 1) such values, the solution Eq. (234) is determined up to a normalization, which 
can be fixed by using Eq. (^281). We thus obtain 
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This equation can be recast for Res < as 



N 



N-l 



(235) 



s'^9h{s) = exp I ^ ln(l - s/pa) - ^ ln(l - 



(236) 



a=l 



The sum over the poles and zeros of the function [0(s) / s^] , with positive real parts, can be rewritten as the following 
complex integrals: 



exp 



{-./. 



-\-ioo 



In 



ln(l - .s/z)| 



4,{z) 



(237) 



This result gives the solution of the problem for an arbitrary bulk Milne kernel. 

The thickness tq of the skin layer can be evaluated by comparing the results Eqs. (235) and (237) with the expansion 
(HI). We get 



To = 



N ^ N-1 ^ 



a=l 



+ioo 



dz 



In 



(238) 



The value for the case of point scattering is obtained by inserting into Eq. ( ^38[ ) the expression for (t>{z) coming from 
Eq. (|225|). By making the change of variable z = itan/3, we get 



To = 



r/2 



d/3 



In- 



tan^/J 



sin^P 3(l-/3cot/3) 



0.710446090 . . 



(239) 



We recover a well-known numerical value (Chandrasekhar, 1960; van de Hulst, 1980; van de Hulst and Stark, 1990). 



B. The inhomogeneous Milne equation 

The above method can be generalized to the inhomogeneous situation (Nieuwenhuizen and Luck, 1993). For point 
scattering one obtains the result 

y TT Jq COs2 /3 + ^2 gin^ (3 J 

This quantity is maximal at normal incidence — 1), where it assumes the somewhat simpler form 



ri(l) = V3 exp j-;^^ ^ df3 ln(l - /3cot/3)| = 5.03647557. . . 
Finally, it is worth noticing that Eq. (|236| ) also implies the following remarkably simple result: 



(241) 



, . n(Ma)n(A^fc) , . 

^(^-^"^^ 3(M. + M.) ■ ^'''^ 

which is particular to the situation in which there are no internal reflections. In the presence of mismatch the 
fmictional form changes. See Eq. (286) for the limit of strong mismatch. 
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C. Enhanced backscatter cone 



Exact results can also be derived for the backscatter cone in the absence of internal reflections. Consider, for the 
sake of simplicity, the enhanced backscatter of a normally incident beam (/ia = 1). The Q dependent Milne equation 
( ^Ol| ) can be solved by means of the Laplace transformation. 

In the case of point scattering, one gets for the backscatter amplitude (Gorodnichev, Dudarev and Rogozkin, 1990) 



1 f 2 r^^ 

7c(Q)=2exp|--y^ < 



d/3 In 1 - 



irctan y/Q^P + tan^ /? 
v/02£2 + tan2 p 



(243) 



For small values of the reduced wave vector Q, the result Eq. ( |243| ) assumes the general form Eq. (221), corresponding 
to the triangular shape of the backscatter cone. The value for Q = reads 



7(1,1) 



= 4.22768104, 



(244) 



in agreement with Eq. (242). In agreement with Table |, one has AQ = l/2£ [see Eq. (221)]. 

We end this section by mentioning the expression for the backscatter cone amplitude for isotropic scattering by 
point scatterers with an arbitrary albedo. 



7c (a; Q) = -exp' 



,^ , / arctan jQ^i^ + tan^ (3 

dp In 1 — a , 

Jo \ y/QH^ + tan2 /3 



(245) 



This result shows that, as soon as the scattering albedo a is smaller than unity the backscatter amplitude is an analytic 
function of Q^, i.e., the cusp at Q = disappears. This confirms the physical intuition that the triangular shape of 
the cone is due to the existence of arbitrarily long diffusive paths. It is therefore a characteristic of the problem with 
unit albedo, i.e. with no absorption, in a half-space geometry. When the reduced wave vector Q and the strength of 
absorption (1 — a) are both small, we observe the following scaling behavior: 



7c (a; Q) 



7c(l;0) 



1 + 2^QWTW 



(246) 



The denominator of this expression is not reproduced quantitatively by the diffusion approximation [see for example 
Eq. (71) of van der Mark, van Albada and Lagendijk (1988)], although it pertains to the long distance physics of the 
problem. We also notice that the prefactor of the square root is nothing but the reciprocal of the value ( |222| ) of iAQ. 
Indeed, Eq. (p4|) yields 7c(0) = 7(1, 1) = t(1)V6, by which Eq. (|22^) reduces to AQ = l/{2£). 



D. Exact solution for internal reflections in diffusive media 

The solution for internal reflections, m ^ 1, can be simply expressed in the solution without internal reflections, 
m — 1. Let us look at T^fP\ the homogeneous solution for the situation with index ratio m. If we define 

then we get from the homogeneous version of Eq. (|152[) 



(247) 



Jo 



We observe that comparison with Eq. (152) immediately gives the solution 

Jo 



(248) 



(249) 



Using the definition Eq. (247) we have to solve the inhomogeneous equation 
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M,) = »^ + ^ / dAW(M,M').H(M'), (250) 
where we used the identities (Nieuwenhuizen and Luck, 1993) 

drrW(r)e--/'' = -Ti(/x) (251) 



/'OC 

Hf^b,l^a)= drT^^^T, 1,^)6-^/'^^. (252) 
Jo 

In a similar fashion we find for the special solution 

/■I 

T'f\T■,^Xa) = r^i\T■,^^a)+ dM5s(M,<^;/^a)4'^(r;M) (253) 

Jo 



with 



9s{n;Ha) = ^7^'H/^;Ma) + ^l^'dn'^^'\n,f,')9s{n';na). (254) 
Following the same steps we obtain the enhanced backscatter intensity. We must solve 

2m V 1 + iQi^Jl - cos 6 / 2/i Jo i-^ 27r 



where (Luck, 1997) 

3(A<a + Mb) 

involves 



(255) 



7«(M«,M.,Q) = ^^^^%^^A^, (256) 



.^^)(,;g)=,^/3exp|-^ ^ ^^^^ , ^ _ arctan y/Q^f^ + tan^\ 1 

' ^ ' \ T^Jo cos2/3 + /z2sin2/? Vg2£2+tan2/32 yj' ^ ^ 

Wc find the following expressions for the dimcnsionlcss injection depth tq, the limit intensity ti, the differential 
reflection dR/dfl, and B{Q), the intensity of the backscatter cone normalized to the diffuse background, at outgoing 
angle Ob coded in Q = k£{6a — 6b)' 



Jrn) _ Jl) 
In 



ri'^ + f diJigH{M^\li) (258) 
Jo 

T^^(Ma) = r['\i,a) + f d/z55(/x;Ma)rf ^(/z) (259) 
Jo 



dOb 47rm2 fiaP-b 



l^"'\p,a,l^b)=l^'^^{lJ-a,IJib)+ d,igs{ii;p,aH'^^{lJi,IJib) (260) 

Jo 
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respectively, 
backscatter. 



The subtracted term in B{Q) comes from single scattering, which does not contribute to enhanced 



For 



Equations (250), (254), and (255) can be solved numerically. In many cases only a few iterations are needed. 
m close to unity the leading corrections for tq and ti(^) are given by the first term in the equation for the g's. 

Explicit results can be obtained in the limit of large index mismatch (Nieuwenhuizen and Luck, 1993). The point is 
that for i?(/x) = 1 the integral kernel becomes {l/2fj,)^^^\fj,, fi'), which has a unit eigenvalue with right eigenfunction 
g{fj.') = 1 and left eigenfunction (7(/i) = 2/i. For a large index mismatch this eigenvalue will bring the leading effects. 
The results will be discussed in the next section. 



E. Exact solution for very anisotropic scattering 



As discussed in Sec. [11 the main effect of anisotropic scattering is that the scattering mean free path i — l/na is 
replaced by the transport mean free path £tr = (cos Q)) — irtr- As a result, the injection depth decomposes as 

zo = Toitr': likewise the opening angle of the backscatter cone is proportional to l/itr- 

For anisotropic scattering the transport equation cannot be solved exactly in general. Several special cases were 
discussed by van de Hulst (1980). He found that in absence of index mismatch the effects of anisotropy are typically 
not very large. Amic et al. (1996) showed that the limit of strongly forward scattering can be solved exactly. This 
extremely anisotropic case is expected to set bounds for more realistic anisotropic scattering kernels. The approach of 
Amic et al. (1996) goes along the lines of previous sections. However, certain constants occurring there are replaced 
by functions of an angle, which are only partly known. For details of the lengthy derivation we refer the reader to the 
original paper. The main results are presented in Table |l[ compared with the case of isotropic scattering. 

We can compare the universal results in the very anisotropic scattering regime, for some of the quantities listed 
in Table ||, with the outcomes of numerical approaches, van de Hulst (1980) has investigated the dependence of 
various quantities on anisotropy for several commonly used phenomenological phase functions, including the Henyey- 
Greenstein phase function 

= n 9 ^~n\ 2W2 - (262) 
(1 — 2gcosB + g )'^' 

The data on the skin layer thickness show that, as a function of anisotropy, tq varies from 0.7104 (isotropic scattering) 
to 0.7150 (moderate anisotropy), passing a minimum of 0.7092 (weak anisotropy); see van de Hulst (1980). The 
trend shown by these data suggests that the universal value 0.718211 is actually an absolute upper bound for tq. 
Numerical data concerning ti(1) is also available, van de Hulst (1997) has extrapolated two series of data concerning 
the Henyey-Greenstein phase function which admits a common limit for very anisotropic scattering {g 1). This 
limit reads in our notation ti(1)/4 = 1.284645, whereas van de Hulst (1997) gives the two slightly different estimates 
1.273±0.002 and 1.274±0.007. The agreement is satisfactory, although one cannot entirely exclude that the observed 
0.8% relative difference is a small but genuine nonuniversality effect. Indeed, the Henyey-Greenstein phase function 
might belong to another universality class than appropriate for the approach of Amic et al. (1996). The same remark 
applies to a less complete set of data (van de Hulst, 1997) concerning the intensity 7(1, 1) of reflected light at normal 
incidence. 

In the presence of index mismatch the very anisotropic scattering limit can be solved exactly along the same lines 
as for isotropic scattering. The formalism becomes very heavy, however. In the limit of a large index mismatch the 
leading behaviors are expected to be the same as for isotropic scattering. The reason is that a large mismatch implies 
that the radiation that has entered the scattering medium will often be reflected internally before it can exit the 
medium. Therefore it undergoes many scatterings, so it loses the nature (anisotropic or not) of the individual events. 
This was verified explicitly by Amic et al. (1996) 



XI. LARGE INDEX MISMATCH 

We have seen that the occurrence of internal reflections complicates the solution of the Schwarzschild-Milne equation. 
Yet in the limit of a large index mismatch the boundaries will act as good mirrors. They re-inject the radiation so 
often that it is also diffusive close to the boundaries. It turns out that, to leading order, the radiative transfer problem 
is greatly simplified. 

Following Nieuwenhuizen and Luck (1993) we consider the Schwarzschild-Milne equation in the regime of large 
index mismatch (m — + , or m ^ 00). In both cases the refiection coefficient i?(/i) is close to unity and the effect 
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of T(/i) = 1 — R{n) is small. We can thus expand in powers of T. We first consider the Green's function of the 
Schwarzschild-Milne equation Eq. (178) and write the kernel as 

M{t, t') = Mb{t - t') + Mi^{t + r') (263) 

= Mb{t - r') + Mb{t + r') - N{t + r') (264) 

N{t + t')= |^'^e-'^T{^,), (265) 



where N is of order J dfiT{ii) ^ 1. This allows an expansion of Eq. ( |178|) in powers of A^. If — (which happens 
for TO — > or m oo), the kernel M has unit eigenvalue, and Eq. (178) has eigenfunction Gs{t) = Cs — const., 
showing that for a medium with perfect mirrors the intensity inside is independent of the conditions outside. For 
small N the value of Cs will depend sensitively on N . Let us write, for iV <C 1, the Green's function as 

Gs(t, t') ^Cs + Go{r, t') + Gi(t, t') (266) 
Cs » 1, Go « 1, Gi « 1. 



(267) 



We insert this in Eq. ( p78| ): 

Gs + Go + Gi = 6{t - r') + Cs 

+ / dr"{AfB(T - r") + Mb(t + t")}(Go + Gi) 
Jo 

poo 

~ dr"7V(r + T")(G5 + Go + Gi). 

^0 

Comparing terms of equal order yields 

Cs = Cs (268) 
Go(T,r') = <5(r-r')+ / dr" {Mb{t ^ t") + Mb{t + t")}Co{t" ,t') 

JQ 

poo 

-Cs / dT"N{T + T") (269) 

^0 

pOO 

Gi(T,r')= / dr"{AfB(T-r")+AfB(T + T")}Gi(T",T') 
Jo 

dr"7V(T + r")Go(T",r'), (270) 



where we have neglected terms of order NGi ~ T^. Integrating these results over r yields an expression for Cs and 
a compatibility condition to ensure that Go has a unique solution: 

Gs = |^ i) ^^'^^"^(^ + ^")} (271) 
drdT"7V(r + T")Go(r", r') = for all r'. (272) 
We can simplify the expression for Cs by introducing the angular averaged flux-transmission coefhcient T : 



OO /'OO 

It 



Jo Jo 

1 N -1 

^d^T(/i) 





Cs^{ I I drdr" I l^e-^T(^) \ (273) 

_ 4 

" r' 



where 



T = 2 fidfiT{fi) = 2 cos 61 sin 61 d6ir(cos 6*). (274) 
Jo Jo 
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Note that /i = cos9 stands for the projection factor of flux in the z— direction. The expression for T can be worked 
out exphcitly using 



m 



-2 



T{^i) = (275) 



For m < 1 we introduce = m ^ — 1 



r = 2/%d.^^lA^^. (276) 



When we substitute 0=sinh (/i/e) this becomes 



, sinh^ (h cosh^ 



-r^&e^l d0 ^ (277) 











2 7o 



2e^ I d(/.sinh2 29!.e-2'^ = ^ / d(j) [e^"^ - 26-2"^ + e"'^"^] (278) 



e-2*+-^— (279) 



Using 



2 2 



sinh(^+ = i = -j=^==, cosh(/)4_ = — ;=^==:. (280) 



TO 



we find exp (2(/)+) = (1 + to)/(1 — to). This yields 

The leading behavior for to ^ could have been found more quickly by taking the leading behavior of the numerator 
and denominator of Eq. (^76| ) . 

rw2/ fidfi^ = -m (to^O). (282) 



Now we do the same for m > 1. We must not forget the Brewster angle: for fi < \J\ — to ^ there occurs total 
reflection. We now introduce = 1 — mr'^ and carry out the integration from e to 1: 

r = 8 /V2dM-^^^2=- (283) 

- ^ ' (to > 1). (284) 



3to2(1 + to)2 

The leading behavior for large to could again have been found by keeping leading terms: 



(m^oo). (285) 



Further we notice that T = 1 for m = 1, as it should be. 

A. Diffuse reflected intensity 

We can now derive simply the leading behavior of the generalized bistatic coefficient: 
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/•OO 

7{fia,f^b)= drdT'G5(r,r')e-^/^''e-^'//'. (286) 
Jo 

dTdr'Cse^^/^^c-^'/^" +C'(T") (287) 



= 4^ ^ ^^^o^_ ^288) 
We can use this to find an approximation for the reflected intensity A^ {da, Ob) in the case m ^ 1. 

cos 9a TaTf, /nan\ 

When m ^ 1, r(/^) ~ Afi/m (since 9 cannot exceed the Brewster angle, which is small at large m), implying 

R cos6'„ 4cos6'a 4cos6'f, 4 ^ 2n a /on^^ 

A — TT ~ cos 9a cos 9b. (290) 

47rm^ m mi -Kin 



B. Limit intensity and injection depth 



The limit intensity Ti(/i) of a semi- infinite space can also be calculated. From Eq. (182) it follows that 



/ dT'Gs(oo,r')e 
Jo 

1 



')e-^'/M 



= lim 



7(M,Aifc) 



r ■ 



(291) 
(292) 
(293) 
(294) 



The injection depth zq — tqI can also be calculated. If we take the (diffusion) approximation Th = tq + r « tq and 
insert it in the above equality we find 



To 



iT' 



Just as in the calculation of Gs we can write the homogeneous solution Th as the sum of three parts: 

rH(T)-^ + ro(r) + ri(T) 

ro(r)=0(l), T^{T)=0{T). 



(295) 
(296) 



(297) 
(298) 



We further require that F p beh aves as r + Too for r ^ oo. The constant Too is a C(l) correction to tq. For m oo 
we can use definition Eq. (265) and that /i w 1 since the Brewster angle is small: 



^ d/i t+t" 
2^ 



r 

e ~T{^l) w — e 



-(r+r') 



7V(t + r') 

Jo 'I 
Inserting this in the equation for Th and comparing terms up to order T^, we obtain 

/■OC 

To = -e-- + / At'[M^{t - t') + Mnir + t')]Fo, 
Jo 



(299) 



(300) 



while the equation for Fi yields the compatibility relation 
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I dre-^ToiT) = 0. 
Jo 



(301) 



These equations cannot be solved analytically. Numerical analysis yields tqq = —1.0357 (Nieuwenhuizen and Luck, 
1993), implying 



exact 
'O 



4 

3T 



1.0357 + 0(r), (m->oo). (302) 



In the limit m — > the transmission coe fficie nt is uniformly small, T(^) sa 3T/i/2. Therefore the term of Eq. (17S) 
arising from the kernel N, defined in Eq. (265), becomes equal to 



X 



d^l'^i' / dTe-^/^'Jo(T,^') = 

d^fi / dTe-^/^ro(T). (303) 



The diffusive shape r + Too is the exact solution of this equation, with tqo = —3/4. This leads to 



4 3 
ST ~ 4 



rr'^' = 1^ - 7 + OiT) (m^O). (304) 



In Fig. |lj numerical values o f ti (1), 7(1 , 1), and 3to have been plotted as a function of m. The solid line is 4/T, 
while the dashed line is Eqs. (|302|) + (^0^. Notice that these asymptotic expressions work quite well up to to = 1. 



C. Comparison with the improved diffusion approximation 



We can compare the exact results for tq with the improved diffusion approximation of Sec. [II. C For m — > 00 we 
can approximate Ci and C2 [see Eq. (|^] as 



Ci = i(l - T) , C2 « i - r MdMr(A^) = 1-1'^- (305) 



From Eq. (|69|) it follows that 



diff _ ^ 
'° ~ 3T 



(306) 



which is very close to the exact value in Eq. ( |302D. No tice that the same result follows when one inserts the diffusion 
approximation Fq (t) = t + Too (for all r) in Eq. ( 301 ) . 

For TO — > deviations from the diffusion approximation show up only to second order in T, 

rdiff = * = ^-l + 0{T) (to ^ 0). (307) 

To find corrections terms in powers of T is essentially as complicated as solving the Milne equation at finite to. For 
TO = 1 (no internal reflections), the exact value is known: Tq^^"^^ = 0.71044. This can be compared with the result for 
the diffusion approximation found in Sec. III.C, 

^cxact ^ 0.71044 (308) 

T^'^ = ^ = 0.666666. (309) 
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D. Backscatter cone 



We now determine the shape of the backscatter cone in the presence of strong internal reflections, in the regime 
m ^ 1. We have to solve 

POO 

Tc{Q,T)^e-^+ / dT'[MB{Q.T-T')+MB{Q,T + T')]Tc{Q,T') 
JO 

poo 

- / dT'N{Q,T + T')rc{Q,T'). (310) 

Jo 

For small Q the solution will have the form 

rc{Q,r)-jc{Q)e''^^- (311) 



Integrating Eq. ( 310 ) over r, using that the bulk-kernel Mb has an eigenvalue and neglecting the Q dependence 

of iV, we find 



yielding 



poo 

= 1+ dTdT'[MB{Q, t-t')- Mb(0, r - r') + 
Mb(Q, t + t')- Mb(0, t + r')]7c(Q)e~'3-' 
- / drdr'iV(0,r + r')7c(Q)e"'5^ 

JO 

f°° 1 T 

= 1 + dr[-l + 1 - -Q2]^c(Q)e-^" - j7c(Q), 

7c(0) = 



r + |g r 1 + Qro q + i/to ' 



(312) 
(313) 

(314) 



(315) 



In passing, we have recovered the triangular shape of the backscatter cone. 

Recently the width of the backscatter cone has been determined experimentally (Den Outer and Lagendijk, 1993). 
The index ratio to was varied using several types of glass containers. The measurements have been compared with the 
theoris of Lagendijk eA, al. (1989), of Zhu et al. (1991) and of Nieuwenhuizen and Luck (1993). The first two results 
have already been discussed in Sec. III.C . The experiments confirmed the theories of Zhu, Pine, and Weitz and of 
Nieuwenhuizen and Luck. 

At first sight it may look strange that an experiment with light (vector waves) can be described by a scalar theory. 
Nieuwenhuizen and Luck (1993) pointed out, however, that for large to only diffusive aspects survive. More subtle 
properties, such as anisotropic scattering and the vector character lead to subdominant effects. The latter was verified 
explicitly by Amic et al. (1997), indicating why the experiments of Den Outer and Lagendijk (1993) are described so 
well by scalar theories. 



XII. SEMIBALLISTIC TRANSPORT 

Wave propagation through waveguides, quantum wires, films, and double-barrier structures with a modest amount 
of disorder can be considered within the same approach. The so-called semiballistic regime occurs when in the 
transverse direction(s) there is almost no scattering, while in the long direction(s) there is so much scattering that 
the transport is diffusive. Various applications were noted by (Nieuwenhuizen, 1993) and discussed in detail by Mosk, 
Nieuwenhuizen and Barnes (1996) The scattering has been considered both in the second order Born approximation 
and beyond that approximation. In the latter case it was found that attractive point scatterers in a cavity always 
have geometric resonances, even for Schrodinger wave scattering. 

In the limit of long samples, the transport equation has been solved analytically including geometries such as: 
waveguides, films, and tunneling geometries such as Fabry-Perot interferometers and double barrier quantum wells. 
The agreement with numerical data and with experiments is quite satisfactory. In particular, the analysis proved that 
the large, gate-voltage independent resonance width of GaAs double-barrier systems, observed by Gueret, Rossel, 
Marclay and Meier (1989) can be traced back to scattering from the intrinsically rough GaAs-AlGaAs interfaces. 
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In the analysis of these systems, the Schwarzschild-Milne equation is decomposed on the cavity modes, thus becoming 
an evolution equation in the z-direction for a matrix. The very same ideas as explained above are used in the solution. 
Due to the complication of the formalism we shall not present details here. We refer the reader to the original works 
and the references therein. 



XIII. IMAGING OF OBJECTS IMMERSED IN OPAQUE MEDIA 

A field of much current interest is the imaging of objects immersed in diffuse media for the purpose of medical 
imaging. The advantage of imaging with diffuse light is that it is non-invasive and as opposed to conventional X-ray 
tomography it does not cause radiation damage. One application is the detection of breast cancer, another application 
is the location of brain regions damaged by stroke (a lot of blood, hence more absorption) or affected by Parkinson's 
disease (less than average blood level, hence more transparency). In these cases one can shine light at many locations, 
and measure the scattered light at many exit points. This information can in principle be used to determine the 
structure of the scattering medium and to detect possibly abnormal behavior. It is therefore important to know the 
scattering properties of idealized objects hidden in opaque media. We shall discuss here some instructive cases. 

For light waves human tissue is a multiple scattering medium like the ones discussed so far. (Yet the absorption 
length is short, typically a few tenths of a millimeter). The diffuse transmission through tissue must already have 
been noticed by our early ancestors when they held their hands in front of a candle. On distances much larger than 
the mean free path, light transport is diffusive and is therefore governed by the Laplace equation. This observation 
brings a deep analogy with electrostatics: an object deep inside an opaque medium is analog to a charge configuration 
at large distance. It is well known that the latter can be expressed in terms of its charge and dipole moment, while 
higher multipole moments are often n(!gligible. A similar situation occurs in an opaque medium: small objects far 
from the surface can be described adequately by their charge (which vanishes if there is no absorption) and dipole 
moment (which describes to leading order scattering from the object). At large distances the effect of, for instance, 
a charge will depend on the geometry of the system. As in electrostatics, mirror charges and mirror dipoles can be 
used in simple geometries. However, to solve the boundary problem of complicated geometries like a head is a hard 
problem. Close to the object, say, within a few mean free paths, the full radiative transfer is needed, making the 
problem more complicated than diffusion theory. 

Den Outer, Nieuwenhuizen and Lagendijk (1993) studied the diffuse image of an objc;c;t cinibeddcd in a homogeneous 
turbid medium and found that it is well described by the diffusion approximation. This framework leads to a prediction 
of the diffuse image of the object, i.e., the profile of diffuse intensity transmitted through a thick slab near the embedded 
object. The shapes of the predicted profiles are universal since they depend only on large distance properties of the 
problem, and were convincingly observed experimentally for the case of a pencil (absorbing cylinder) and a glass fiber 
(transparent cylinder). Further results on imaging of objects within the framework of the diffusion approximation can 
be found in Schotland and Haselgrove (1993), Feng, Zeng and Chance (1995), Zhu, Wei, Feng, and Chance (1996), 
and Paasschens and 't Hooft (1998). 



A. Spheres 

In the diffusion approximation the problem is straightforward. Close to the object (though still many mean free 
paths away) the intensity can be written as 

/(r) = /o(r) + <57(r), (316) 

where /o(r) is the intensity profile in the absence of the object. Luck and Nieuwenhuizen (1998) have introduced the 
following notation for the disturbance of the intensity 

57(r) = (g + p.V)^— = -p--^^^. (317) 

|r-ro| |r-ro| |r-ro|'^ 

q is called the "charge" and p the "dipole moment" of the object. The linearity of the problem and the isotropy of 
the spherical defect imply that the charge q is proportional to the local intensity /o(ro) in the absence of the object, 
while the dipole moment p is proportional to its gradient. We therefore set 

q = -Q/o(ro), P = -PV/o(ro), (318) 
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Q is called the capacitance of the spherical object, and P its polarizability. These two parameters are intrinsic 
characteristics of the object. The capacitance Q is non-zero (and positive) only if the sphere is absorbing. 

For plane wave incidence on a slab geometry one has the familiar expression /o(r) — {L — z)Iq/L. A double set of 
mirror charges is needed to make I — {L — z)Iq/L vanish at z = and z ^ L for all p: 



T oo , 

n— — oo ^ 



1 1 



[(z - zo + 2nLY + {p-paY] [{z + zo + 2nL)2 + [p - p^f] ^1' 



(319) 



El z - zq + 2nL z + zq + 2nL i 
1 [(z-zo + 2nL)2 + (p-po)2]3/2 + [(z + zo + 2nL)2 + (p-po)2]3/2 1 • ^ ' 



n— — oo 



The intensity T{x,y) transmitted through the sample, and emitted on the right side (z = L) at the point p ~ {x,y), 
is proportional to the normal derivative of the diffuse intensity at that point: 



T[x,y) = -Ki 



(321) 



dz 

When there is no embedded object, the transmission is uniform To = KIq£/L. In the presence of an object one finds 

L — zo + 2nL 



r(p)-ro i-2g(i-zo) J2 



+ 00 



n=-oc. [{L - Zo + 2nL)2 + (p - po) 
2{L -Z0 + 2nLf - {p - p^f 



3/2 



-2P y -l---o^-"-^ ■ (322) 

nt^oo [{L - ZO + 2nLY + {p- poY] ) 

The information on the type of scatterer is coded in its capacitance Q and polarizability P. In experiments on pencils 
and fibers, 2 — d analogs of these profiles have been clearly observed experimentally ( Den Outer et ai, 1993). In the 
same fashion the reflected intensity can be estimated by taking the derivative in Eq. ( |32l[ ) at z = 0. The rest of this 
section is devoted to determining Q and P in idealized situations. This amounts to solving the scattering problem 
with the radiative transfer equation in the presence of an object placed at ro. 

For large objects the diffusion theory can be used which simplifies the problem. Following Den Outer et al. (1993) 
we consider a macroscopic object of radius R, having diffusion coefficient D2 and inverse absorption length k, while 
the medium is nonabsorbing and has diffusion coefficient Di. At the boundary of the object one has 



dn 

Near the object one has 



Iout{R+)^hn{R-) (323) 

(324) 



„ diout 



D2 



R+ On 



z z 

lin = Asinhc(Kr) :^[sin\ic{nr) — cosh(Kr)] , (325) 



where sinhc(a;)=sinh(a;)/x. Matching the solutions gives 



Q ^ ^ £'2 (cosh kR - sinhc kR) 

Disinhc kR + D2 (cosh kR — sinhc kR) 



For a totally absorbing sphere (k 00) the result Q = R follows immediately by requiring in Eq. ( |316| ) that I{R) = 0. 
The polarizability reads 

p ^3 Z?i (sinhc Ki? — cosh Ki?) — £'2(2 cosh Ki? — 2sinhcK_R — Ki?sinhKi?) (327) 
2i?i (sinhc kR — cosh kR) + D2 (2 cosh kR — 2sinhc kR — kR sinh kR) 

In the absence of absorption this reduces to 



,3 Di - D2 



Q = 0, P-R^^j^^. (328) 
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Lancaster and Nieuwenhuizen (1998) consider the limit of small objects: a point scatterer with scattering cross 
section age, extinction cross section aex^ and albedo Oe = Csc/cex- If diffuse intensity hits the scatterer once, but does 
not return to it, Q and P read 

Q^ = l'-f. (329) 

where aats — cTex — <^sc- The intensity can also return an arbitrary number of times to the scattererer. For point 
scatterers this boils down to a geometric series expansion for Q and P. 

where Bi and B2 are integrals describing the propagation of the diffuse intensity through the medium and the effect 
of hitting the scatterer another time. These integrals are actually divergent in the limit k£ ^ 00, which leads to 
introduce extra internal parameters of the point scatterer. 

As opposed to extended objects, for point scatterers the maximally crossed diagrams also contribute. The reason 
is that exit and return point coincide, so there is no de-phasing. The analysis finally yields 

Q = oT' P= oT' (331) 

1 + Qi (2Bi - b[^> ) 1 + Pi (2^2 - P^'^ ) 



which retains the structure of Eq. ( |330D . Note that there is an additional factor 2, except for the diagrams with one 
intermediate common scatterer, an effect that is accounted for by subtracting b[^1. This is very similar to the fact 



that the single-scatterer diagram does not contribute to the enhanced back scatter cone, see eq. ( p09[) . Lancaster and 
Nieuwenhuizen (1998) also considered more realistic situations in an approach involving a partial wave expansion. 

For extended objects it is a highly non-trivial task to go beyond the diffusion approximation. Luck and Nieuwen- 
huizen (1998) considered the application of the radiative transfer equation to such cases. One thus has to solve the 
analog of the scattering problem of quantum mechanics, for this transfer equation. This is generally a hard and still 
unsolved problem. However, there are three cases where the inside of the scatterer plays a trivial role: a totally 
absorbing (black) object (having k = 00), a transparent, non-absorbing object (having k — 0, D2 = 00), and a totally 
reflecting, non-absorbing object (having k = D2 = 0). For spheres and cylinders of these types, the scattering problem 
can be reduced to an integral equation in one variable, which is solved exactly in the limits of small and large objects. 
In the latter case the results of diffusion theory are recovered, and the first order corrections to diffusion theory were 



found, see Table III and Figs. 13 and 14, The same approach can be extended, in principle, to more general cases. 



B. Cylinders 

The use of the radiative transfer approach can be extended to other objects (Luck and Nieuwenhuizen, 1998). 
Consider a cylinder with radius R with its axis parallel to the x-axis. Assume the sample is infinitely long in the x- 
direction, so that both Io(r) and the disturbance 51{r) depend only on the two-dimensional perpendicular component 



= (y, z). Eq. (317) is replaced by 



«(r^) = (g-fp.V)ln^^ = (7ln— (332) 
|r-ro| F-ro| |r-ro|^ 

The length scale R, requested by dimensional analysis, is determined by conditions at the boundaries of the sample, so 
that it is in general proportional to the sample size. This sensitivity of SI{r) to global properties of the sample arises 
because the logarithmic potential of a point charge in two-dimensional electrostatics is divergent at long distances. 
As a consequence, the capacitance Q is not intrinsic to the cylinder. An intrinsic quantity is its effective radius 
RcS, defined by the condition that the cylindrically symmetric part of the total intensity, including the charge term, 
vanishes at a distance r — R^s from the axis of the cylinder. We thus have /o(ro) + 9 ln(-R/i?eft) — 0, hence, with 
q = -Ii-ro)Q, 

g=^^. (333) 
In^ 

Half 
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We observe that ReS and R have the dimension of a length, while the polarizability P has the dimension of an area. 
As a consequence, the orders of magnitude -Roff ~ R and P ^ R^ can be expected for a cylinder of radius R. For 



the same physical situations as above, more accurate predictions can be found in Table [II. For a cylinder in a slab 
geometry the transmitted intensity T{y) finally reads in terms of Q and P 

( / zo\ sin^SL ^2 i + cos^cosh^\ , ^ 

\ ^ \ cos + cosh ^ ^cos2£o +cosh^)V 



where 



R~ — sm (335) 

TT L 

is indeed of order L, but also depends non-trivially on the aspect ratio zq/ L. These shapes have been observed in the 
experiments of Den Outer et al. (1993). 



XIV. INTERFERENCE OF DIFFUSONS: HIKAMI VERTICES 

So far our discussion has focused on the intensity aspects of wave transport. We now consider interference effects, 
which are due to the inherent wave nature of the problem. First we have to introduce the concept of interference 
vertices. These vertices describe the interaction between diffusons. The vertices are named after Hikami, who used 
them in 1981 (Hikami, 1981) for the calculation of weak-localization corrections to the conductance. Their original 
introduction dates back two years earlier to the work of Gor'kov, Larkin and Khmel'nitskii (1979). Remember that 
each diffuson consists of two amplitudes. The vertices describe the exchange of amplitudes between diffusons, which 
leads to correlations between the diffusons. It is also possible to represent these processes using the standard impurity 
technique diagrams, as described for instance by Abrikosov, Gor'kov and Dzyaloshinski (1963), but the interference 
vertices prove to be more convenient. 



A. Calculation of the Hikami four-point vertex 

We start with the simplest vertex: the four point vertex, in which two diffusons interchange an amplitude. The 
technique for calculating the vertex by means of diagrams is well known (Abrikosov et al., 1963): First, draw the 
diagrams; second, write down a momentum for each line and use momentum conservation at each vertex; and finally, 
integrate over the free momenta. The first step, writing down all leading diagrams, requires some care. In Fig. ^ we 
have drawn the diagrams in the Hikami box to second-order Born approximation. It is important to include all leading 
diagrams, as there will be a cancellation of leading terms. These cancellations are imposed by energy conservation. In 
the second-order Born approximation one neglects scattering more than twice from the same scatterer; the (dashed) 
interaction line indicates that the propagators scatter once on the same scatterer. Scattering processes involving just 
one propagator were already included, as the propagators are dressed. We have not drawn all diagrams, one can 
imagine. First, the box is attached to a diffuson ending with a scatterer (which is consistent with the way we have 
defined the diffusons). As a result, common scatterings between propagators on the same leg are not allowed in the 
vertex as they are already included. Secondly, diagrams with two or more scatterers and parallel dashed lines are 
subleading. They contain loops with two propagators of the same type, i.e., integrals like J d'^pG'(p -I- q)G(p), which 
turn out to be of higher order in l/(fco^)- Third diagrams with two crossed dashed lines are also subleading. Finally, 
we again work in the independent-scatterer approximation. 

As an example we calculate the second right hand side diagram of Fig. |l^. As usual in the diffusion approximation, 
the diffuse intensity is assumed to be slowly varying, which allows us to expand the Green's functions as 



+ 2p ■ q + - — - Alo— - nt 



G - (2p • q + _ Ac^-)G2 + 4(p • q)2G3, (336) 



with the shorthand notation G = G(p,a;). The momenta, numbered according to Fig. |l^, point towards the vertex. 
We number the Green's functions also according to the figure. Note that this particular diagram is a product of two 
loops. To lowest order in we have 
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C 



X [h 



(337) 



where the q^ are the momenta of the diffusons, the oji are the frequency differences of the diffusons, and u is the 
scatter potential. The /— integrals are given in Appendix |A[ We take absorption terms into account to lowest order, 
i.e., only in the leading contributions, the /2,i terms. This yields 
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(-q2 • q3 - qi • q4 + X! + + 



-it 



(Wi + W3 — ^2 — W4). 



(338) 



Here we made use of the fact that in second order Born approximation nii} — 47r/£. 

A similar calculation gives n'f^ and Defining the reduced frequencies of the legs as J^i = — 104) /D, 

0.2 — — (wi — L02)/D, r23 = — (cJ3 — L02)/D and Vt^ = — {uj^ — ijj4)lD, we find for the sum of the three diagrams 



Hi = H 
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i +^4 
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-2qi-q3-2q2-q4 + ^(q- 



iVLj 



(339) 



This is the main result of this section. Note that the leading, constant terms proportional to (? jk'^ have canceled, 
which is closely related to the current conservation laws (Kane, Serota and Lee, 1988). 

In Eq. ( |339| ) it is important to keep track of the terms. When the vertex is attached to a diffuson, the (j| terms 
together with k and 17 yield according to the diffusion equation, (q^ ^ ^ inj)Cj ~ 12Ti/t. The constant factor in 
the right-hand side corresponds to a delta function in real space. If we attach external diffusons to the Hikami box, 
this delta function is nonzero at roughly one mean free path from the surface and its contribution can be neglected, 
as its effect is of order ^/ L. As a result, the absorption part and the frequency-dependent part vanish: 



-ff4(qi,q2,q3,q4) = -/i4(qi- qa + q2- q4)'5(Eq«) ' 



(340) 



where we have denoted /14 — t / (487rfc2). In most calculations, however, we use a Fourier transform in the z direction, 
because in our slab geometry the {qx, qy, z) — (qj_, z) representation is the most convenient. Then the qz terms become 
differentiations: 
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(341) 



The differentiations work on the corresponding diffusons and afterwards Zi should be set z. Formally one has 

Hi{z; Zi,Z2, Z3, Zi)Ci{zi)C2{z2)C3{z3)C4{z4) 

hi 



[2dz,dz 



Ci{zi)C2{z2)C3{z3)Ci{zi)\^ 



(342) 



but we just write this as Hi{z)£i(z)C2{z)C3{z)Ci(z). Note that the vertices yield the spatial derivatives of the 
diffusons, that is, their fluxes. This observation is helpful in estimating the influence of internal reflections, which 
reduce the spatial derivatives of the diffusons. 



B. Six-point vertex: 

We shall also need a second order diagram the six-point vertex Hq. Six diffusons are connected to this diagram 
(Fig. |l^). This diagram was calculated by Hikami (1981). Again, the dressings of the diagrams have to be added to 
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the bare diagrams. Taking rotations of the depicted diagrams into account, there are 16 diagrams in the second-order 
Born approximation. (It is not allowed to dress the bare six-point vertex with a scatterer that connects two opposite 
propagators. This dressing gives a leading contribution also, but the resulting diagram is the same as a composed 
diagram with two four-point vertices. Therefore, it should not be included here, but it will enter in diagrams with 
two H4 vertices.) One finds 



967rfc4 

6 



[qiq2 + q2q3 + q3q4 + q4q5 + q5q6 + q6qi 
1 



2 



(343) 



Apart from the frequency and absorption terms included here, Hikami's original expression can be derived from 



Eq. (343) using momentum conservation. 



C. Beyond the second-order Born approximation 

Going beyond the second-order Born approximation in the calculation of the diffuson diagrams, meant that in the 
diagrams the t matrix replaced the potential u. This resulted in a replacement of the mean free path: it became 
£ = nttlA'K instead of £ = nu^ /Att. However, for the vertices the situation is more subtle. In the previous calculations 
we worked to second-order of the scattering potential, thus neglecting diagrams with more than two scatterings on 
the same scatterer. Including all orders there are eight diagrams. Fig. ^ 

The calculation is very similar to that of the second-order Born approximation above. The only difference is that 
instead of m, the t matrices t and t occur {t on G"s, t on G*'s). After summing the expressions for the eight diagrams. 



one finds the same result for the Hikami box as was found previously in second order Born approximation, Eq. (339). 
However, the definition of the mean free path is different: I = ntt/An instead of ^ = nv? /A-k. We conclude that, 
although extra diagrams were present, only the mean free path was renormalized in extending from second-order Born 
to the full Born series. This can also be shown on a much more generally in the nonlinear sigma model (Lerner, 1996). 
For the six-point vertex it was not checked explicitly whether the second and full Born approximations give the same 
result (there are at least 64 diagrams), but we have no doubt they will. 



D. Corrections to the conductivity 

In the original works of Gor'kov et al. and Hikami, the vertices were introduced to calculate weak-localization 
corrections to the conductivity. We have already mentioned that the maximally crossed diagrams give the largest 
correction to the conductivity, increasing the return probability of the intensity. In the Hikami vertex formalism these 
diagrams are drawn in Fig. |l8|. The correction to the diffusion constant diverges with the sample size in one and two 
dimensions. This divergence indicates the absence of extended states in one and two dimensions no matter how small 
the disorder. In three dimensions, however, there is a transition to a localized state at a finite level of disorder. In 
this work we do not consider such loop effects. We suppose that we are far from localization, so that we can restrict 
ourselves to the leading processes, and the vertices show up only in the interaction between different diffusons. 



XV. SHORT RANGE, LONG RANGE, AND CONDUCTANCE CORRELATIONS: Ci, C2, AND C3 

A nice feature of optical mesoscopic systems, as compared to electronic systems, is that they permit to measure 
certain parts of the transmission signal separately. Depending on whether integration over incoming or outgoing 
directions is performed, one measures fundamentally different quantities. In electronic systems one usually measures 
the conductance of the sample. This is done by connecting incoming and outgoing sides to a clean electron "bath." 
All electrons are collected and all angular dependence is lost. In contrast, in optical systems one usually measures 
the angular resolved transmission, but angular integration is possible. The correlation functions of the different 
transmission quantities, their magnitudes, their decay rates, and the underlying diagrams are all very different. 

In this section we discuss correlation functions for all three transmission quantities. See also the review by Berkovits 
and Feng (1994) for the physical meaning of the correlations. We try to provide all quantitative factors relevant for 
an experiment. In later sections the same approach will be extended to higher order correlations. 
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A. Angular resolved transmission: the speckle pattern 



If a monochromatic plane wave shines through a disordered sample, one sees large intensity fluctuations in the 
transmitted beam, so-called "speckles." The speckle pattern is wildly fluctuating (as a function of the frequency of 
the light or the outgoing angle). The typical diameter of a speckle spot on the outgoing surface is A. We are interested 
in the correlation between speckle patterns of two different beams. The beams may have different incoming angles, 
different frequencies, or different positions. 

We denote the transmission from the incoming channel a (wave coming in under angles 9a, 4>a) to the outgoing 
channel 6 (waves transmitted into angles Oh, 4>b) as Tab- There are short, long, and "infinite" range contributions to 
the correlation function of the speckles. The frequency correlation functions can be classified as 

{Tabiu;)T,,{u; + Au;)) ^ ^ ^ cf'd^^^^ ^ Cf-^{Au;) + Cf^^iAuj). (344) 



{Tab{L0)){T,di^ + AL0)) 



Feng, Kane, Lee and Stone (1988) originally put this expression forward as an expansion in 1/g, where Ci was of 
order one, C2 of order 1/g, and C3 of order We return to this expansion below. We present a sketchy picture of 

the diagrams of the different correlations in Fig. |l^. The largest contributions to transport are from the independent 
diffusons, yet correlations are present, mixing diffusons with different frequencies or angles. 



The unit contribution in Eq. (344) just comes from the uncorrelated product of average intensities. The Ci-term in 
the correlation function is the leading correlation if the angular resolved transmission Tab is measured. It is unity if 
the two incoming directions a and c are the same, the outgoing directions b and d are the same, and Alu is zero. If one 
changes the frequency of the incoming beam, the speckles of the outgoing beam deform and eventually the correlation 
vanishes exponentially. This effect is the short ranged Ci contribution. Also of interest is the case of the angular Ci 
correlation with the frequency fixed: If the angle of the incoming beam is changed, the speckles change due to two 
effects. First, the speckle follows the incoming beam. This is also known as the memory effect, which one can even 
see by the naked eye. Secondly, the speckle pattern also deforms, i.e. it de-correlates, this is again the Ci correlation. 
The Ci correlation is a sharply peaked function nonzero only if the angle of the incoming and outgoing channel are 
changed by the same amount. Diagrammatically, the Ci factorizes in two disconnected diagrams, see Fig. 

Theoretical studies of Ci were first done by Shapiro (1986). Experimental work on Ci as function of angle, explicitly 
showing the memory effect, was done by Freund, Rosenbluh and Feng (1988). van Albada, de Boer and Lagendijk 
(1990) studied the frequency dependence of Ci. Later, the effects of absorption and, especially, internal refiection 
were studied (Freund and Berkovits, 1990; van Rossum and Nieuwenhuizen, 1993). 



B. Total transmission correlation: C2 



The C2 correlation stands for the long range correlations in the speckle pattern. Consider again a single-direction-in 
single-direction-out experiment. Again the two incoming angles are the same {a — c), yet this time we look at the 
cross correlation of two speckles far apart. As there is an angle difference in the outgoing channel, but not in the 
incoming beam, the Ci term is now absent. Instead there is a much weaker correlation. As the frequency shifts, this 
correlation decays algebraically. This is the C2 contribution. It describes correlations between speckles far apart. 
In a single-channel-in single-channel-out experiment it is only possible to see the weak effects of these higher order 
correlations in very strongly scattering media, i.e., if the system is rather close to Anderson localization. This type of 
experiment was done by Genack and Garcia (1993). The C2 term can be measured more easily in a setup using one 
incoming channel and integrating over all outgoing channels. In an experiment the outgoing light is collected with 
an integrating sphere. Therefore, one measures the total transmission Ta = X^b^afe- this setup the C2 correlation 
function, which has a much smaller peak value but is long ranged, contributes for all outgoing angles. The C2 term is 
now the leading correlation as the sharply peaked Ci term is overwhelmed, a precise estimate of the Ci contribution in 
C2 was given in (van Rossum, de Boer and Th. M. Nieuwenhuizen, 1995). Only outgoing diffusons with no transverse 
momentum and no frequency difference in their amplitudes are leading in the total transmission, as then the phases 
of the outgoing amplitudes cancel. From Fig. one sees indeed that it holds for the C2 diagram (outgoing lines of 
similar style pair), but not for the Ci diagram. The long range character of the C2 arises due to interference of the 
diffuse light paths. It is of order g^^ and corresponds to a diagram where the two incoming diffusons interact through 
a Hikami vertex, where the diffusons exchange amplitudes. 

The long range C2 correlation function was first studied by Stephen and Cwilich (1987). Zyuzin and Spivak 
introduced a Langcvin approach to simplify the calculation of correlation functions (Zyuzin and Spivak, 1987). Pnini 
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and Shapiro (1989) applied this method to calculate the correlation functions of light transmitted through and reflected 
from disordered samples. The C2 correlation functions were measured in several experiments. For optical systems 
Van Albada et al. performed measurements (van Albada et ai, 1990; de Boer, van Albada and Lagendijk, 1992). 
Microwave experiments were performed by Genack, Garcia and Polkosnik (1990) and Garcia, Genack, Pnini and 
Shapiro (1993). In several papers effects of absorption and internal reflections were studied, see Pnini and Shapiro 
(1991), Lisyansky and Livdan (1992), Zhu et al. (1991), and van Rossum and Nieuwenhuizen (1993). It was shown 
that absorption and internal reflection, neglected in the earliest calculations, signiflcantly reduce the correlations. 



C. Conductance correlation: C3 



Finally, the C3 term in Eq. (344) is dominant when the incoming beam is diffuse, and one collects all outgoing light. 
Experimentally one could do this using two integrating spheres. Then one measures, just as in electronic systems, 
the conductance g = (T) = b^ab- In that measurement contributions to the correlation are dominant, where the 
angles a and c are arbitrary far apart. Though C3 is of order g^^, it dominates over the Ci and C2 terms as it has 
contributions from all incoming and outgoing angles. Therefore it is sometimes called the "infinite" range correlation. 
In contrast to the C2 term, now also the amplitudes of the incoming diffusons must have opposite phases. This occurs 
in a diagram where the two incoming diffusons interact twice. Note that a loop occurs, see Fig. ^| Apart from the 
normalization to the average, C3 equals the Universal Conductance Fluctuations, or UCF. 

The UCF were an important discovery in the study of mesoscopic electron systems. The electronic conductance of 
mesoscopic samples shows reproducible sample-to-sample fluctuations. This was observed in experiments on meso- 
scopic electronic samples by Umbach, Washburn, Laibowitz and Webb (1984). The theoretical explanation was given 
by Altshuler (1985) and by Lee and Stone (1985). For reviews on the subject see Lee, Stone and Fukuyama (1987) 
and the book edited by Altshuler et al. (Altshuler, Lee and Webb, 1991). The discovery of these fluctuations showed 
that interference effects are important in electronic systems, even far from the localization transition. The variance 
of the fluctuations is independent of the sample parameters such as the mean free path, the sample thickness and, 
most remarkably, the average conductance. Therefore, the fluctuations are called universal conductance fluctuations. 
The conductance fluctuates when the phases of the waves in the dominant paths change. This happens, of course, if 
one changes the position of the scatterers, e.g. by taking another sample. One also may keep the scatterers fixed but 
apply a magnetic field or vary the Fermi energy. In all these cases one modifies the phases of scattered waves, so that 
different propagation paths become dominant. The fluctuations are much larger than one would obtain classically 
by modeling the system by a random resistor array, in which interference effects are neglected. A classical approach 
is valid only on a length scale exceeding the phase coherence length, where the fluctuations reduce to their classical 
value. 

Unfortunately, despite some advantages of optical systems, the analog of the UCF has not yet been observed in 
optical systems. Such experiments turn out to be difficult. Although the magnitude of the fluctuations is universal, 
they occur on a background of order g, where g is the dimensionless conductance (in optical experiments one typically 
has g ~ 10^). The relative value of the fluctuations to the background is thus 1/g, so that the C3 correlation function 
is of order l/.g^, typically of order 10^^. For electrons this problem is absent as smaller values of g are achievable. In 
the electronic case there have been many observations of the universal conductance fluctuations. 

A promising candidate for measuring C3 is the microwave experiment of Stoytchev and Genack (1997). In such a 
set up one always has good control over the channels. Moreover, one is able to study materials with large disorder, 
i.e. down to g = 3. 

In principle the C2 and C3 correlations are present as a sub-leading term in the angular correlation function, 
but there are also other contributions of order 1/g and l/g^, respectively, to the correlation Eq.(344). The weak 



localization correction, see Sec. XIV. D, is an example of a 1/5 contribution to Ci. We have drawn a 1/g contribution 
to C2 in Fig. Such corrections are studied by Garcia et al. (1993). Therefore, the interpretation of the C's as 
an expansion in 1/g is misleading. Rather we define the Ci, C2, and C3 as the leading term in the correlation of 
Tab, Ta, and T, respectively. Equivalently, Ci contains disconnected diagrams; incoming and outgoing amplitudes 
automatically have the same pairing. The C2 term contains connected diagrams swapping the initial pairing, and C3 
contains all connected diagrams in which incoming and outgoing pairings are identical. 



XVI. CALCULATION OF CORRELATION FUNCTIONS 
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A. Summary of diffuse intensities 



For the calculation of interference effects it is often too complicated to take the precise behavior near the boundaries 
into account. We therefore approximate the diffuse intensity with tis bulk behavior, yet we will use the prefactors 
and angular dependence, i.e. especially tq and ri, found with the radiative transfer approach. 

The diffuse intensity created from a beam impinging in direction a is 



k£A^a L 

For uniform illumination of the sample the source term in the transport equation becomes 



(345) 
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Using this source term in the Schwarzschild-Milne equation gives the diffuse intensity for a uniform illuminated sample 

Ak L + zq — z 
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Defining 



L + 2zQ 
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and comparing the pre-factors, one has 

Finally, for diffuse propagation between two points in the slab one has 

—^Cint{z,z';M) + M^Cintiz,z';M) = ^5{z-z'), 



with the solution 



Cint{z,z';M) = 



127r sinh M z< sinh M {L - Zy) 



M sinh ML 



For angular resolved measurements the outgoing intensity in a certain direction b is, in analogy with Eq. (349), 



The angular transmission coefficient reads 



(Tab) 

The total transmission coefficient is 
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The conductance reads, cf. Eq. ( |190| ) 



(T) = J d^rCout{z)S{z) 
9kA t°° 
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(356) 
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B. Calculation of the Ci correlation 



The calculation of Ci is rather simple as it is just the product of two independent diffusons with exchanged partners. 
Both diffusons consist of one amplitude from one beam and one complex conjugated amplitude from the other beam. 
Due to momentum conservation in each diffuson, the perpendicular momentum difference of the two incoming beams, 
q^, is equal to the difference of the outgoing momenta. Therefore, the Ci is only nonzero if incoming and outgoing 
angles are shifted equally during the experiment. 



We calculate the normalized product of a diffuson with "squared mass" 
conjugated diffuson. This yields: 



iil and the complex 
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k'^Zq) sinh kL 
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In Fig. ^this function is drawn for various situations. In thick samples where Mzq ^ 1 and no absorption is present, 
the result reduces to (Feng et ai, 1988) 



Ci(M) 



ML 



sinh AIL 



(358) 



For large angle or frequency difference (|M|L 1), Ci decays exponentially, as Ci(M) ~ exp(— 2Re(M)L). 
spatial correlations of beams can be extracted from the Fourier transform of the angular correlation function. 



The 



C. Calculation of the C2 correlation 

The diagram of the long range correlation C'2 is depicted in Fig. Remember that Ci is zero if incoming and 
outgoing angles are changed unequally, but C2 is connected, allowing momentum flow from one diffuson to the 
other, therefore it is nonzero if one changes incoming and outgoing angles by a different amount. C2 stems from the 
interaction between two diffusons, which exchange partners somewhere inside the slab. The Hikami box describes 
this exchange. (For an estimate of Ci correlations a C2-type experiment see van Rossum, de Boer and Th. M. 
Nieuwenhuizen (1995).) We have 



1 



L 



C2 = / dzCUz, Mi)^^{z, M3)Hi{z)Cout{z, M2)£out(^, M4) , (359) 

{-l-a){J-c) Jo 

We have labeled the incoming diffusons with odd numbers and the outgoing diffusons with even numbers, a convention 
used throughout this work. Attaching diffusons to the Hikami box allows a simplification of the expression for the 
box. We calculate C2 by inserting MI4 = q±'^ + ± Mi = M3 = k. In the z— coordinate representation the box 



reads, see Sec. XIV 



H4{zi,z2,z3, Z4) = h4{dz-,dz3 + dz^dz^ + q±'^) . (360) 
The above calculation for C2 is similar to the one obtai ned i n the Langevin approach, as one can see by comparing 



Eq. (39) of (Pnini and Shapiro, 1989) to Eqs. (|36C| ) and ( 359 ) in this section. In the Langevin approach one assumes 



that there is a macroscopic intensity, which describes the average diffusion, and an uncorrelated random noise current 
superposed. In our approach we have a slowly varying diffuse intensity, while the point-like, uncorrelated random 
interactions originate from the Hikami box. 

We discuss some simplified cases. Consider first the case of angular correlations and no absorption but taking the 
boundary effects into account: 

C2{q^)^j^F2{q^L), (361) 
in which we define the dimensionless function F2 
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F2iq±L) = [smh2q^L - 2q^L + q^za (6 sinh^ q^i - 2g±^L^) 

6q±^ZQ sinh^ q±L - 
[2qj_L {{1 + q±^ zl) smhq±L + 2q±ZQ coshq^lYy 



V^qi_^ZQ (sinh2(7^L — q^L) 



sinh2g^i] X 



(362) 



which decays like \/q± for large q± {i.e. for large angles). However, it still should hold that q^^^ <^ 1, else we 
probe essentially interference processes within one mean free path from the surface, which is beyond this approach. 
Neglecting boundary effects, zq <C i, one recovers (Pnini and Shapiro, 1991) 



F2{q^L) 



sinh(2g^i) — 2q^L 
2q±L sinh^(gj^L) 



(363) 



with ^2(0) = ^. In Fig. ^ we have plotted the F2 function for various situations. 

We have studied the correlation as function of the two-dimensional momentum q± . Similarly, we can study the real 
space correlations. The functional form in real space is of roughly similar shape as Fig. |2^, with typical decay length 
L in the x, y direction (Stephen and Cwilich, 1987). One can also obtain frequency correlations. Without absorption, 
neglecting boundary effects, one finds (de Boer et ai, 1992) 



2 I sinh( V 2fi L) - sin( V 2fi L) 
Vcosh(y2^-/^) - cos(y2^L) 



(364) 



Of special interest is the top of the correlation function. By definition the top of the correlation function corresponds 
to the second cumulant, or second central moment, of the total transmission distribution function. With a plane wave 
as incoming beam, all transverse momenta are absent, we also neglect absorption. These assumptions correspond 
to all M being equal zero in Eq. ( |359|) . The diffusons are simple linear functions, given by Eq. (345) and ( 352 ). 
Neglecting internal reflections, one obtains for the second cumulant (Stephen and Cwilich, 1987) 
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where double brackets denote cumulants normalized to the average. The correlation is thus proportional to I/5. The 
inverse power of g counts the number of four-point vertices, this power can be looked upon as the chance that two 
intensities interfere. 



One obtains a general expression for C2 by inserting the appropriate diffusons and the Hikami box Eq. (360) into 
Eq. (|359|). The result is rather lengthy and was given in van Rossum and Nieuwenhuizen (1993). 



1. Influence of incoming beam profile 

We now study the influence of the beam profile on the correlation. In experiments the beam is often focused 
to a small spot in order to minimize the dimensionless conductance g and therefore to maximize the correlations. 
Physically it is clear that the correlations increase if the two incoming channels are closer to each other, i.e., if the 
beam diameter is smaller. If the spot of the incoming beam is finite, amplitudes with different angles, i.e. different 
transverse momenta, are present. They can combine into incoming diffusons with perpendicular momentum. We 



suppose that the incoming beam has a Gaussian profile. We decompose it into plane waves defined in Eq. (147) (for 
convenience we assume the average incidence perpendicular) 

= § E -^(^^aX . 0(9i) = -|=e-'-^''«/^ (366) 

where po is the beam diameter. In order to have two diffusons with a momentum q±^ and — g^, we find that the four 
incoming amplitudes combine to a weight function 



56 



J d^Fid^Pg 0(Pi)0*(Fi + q±) cj){P3)cj)*{P3 - q^) = ^-pIi^"I\ (367) 

We get the correlation function by integrating the dependent correlation function over the momentum with the 
corresponding weight. Neglecting boundary layers we find for the top of the correlation (de Boer et ai, 1992) 



(368) 



with F2 as in Eq. (363). If the incoming beam is very broad, po ^ L, only the term F2{q± = 0) = 2/3 contributes 
and one recovers the plane wave behavior {{T^)) — 2/3g. The second cumulant decreases as 1/pq at large po. At 
smaller beam diameters the top of the correlation is proportional to 1/poi ^ind diverge in this approximation, in which 
q^e - £/po < 1. 



2. Reflection correlations 

Similar to transmission correlations, reflection correlations were studied. We have seen that the transmission corre- 
sponds to the diffusons, but in reflection there are two intensity diagrams: the diffusons bring a constant background, 
the maximally crossed diagrams yield the enhanced backscatter cone, which is only of importance if incoming and 
outgoing angles are the same. This leads to more diagrams and to more peaks in the reflection correlation than in the 
transmission correlation. For the Ci one can pair diffusons and maximally crossed diagrams (Berkovits and Kaveh, 
1990). Indeed in optical experiments two peaks were seen by Freund and Rosenbluh (1991). The C2 term in reflection 
was treated by Berkovits (1990). Measurement of C2 in reflection by collecting all reflected light is tricky as it is hard 
not to interfere with the incoming beam. Nevertheless, measurement may be possible as a long range component in 
the angular resolved reflection (this however requires quite small values of g). For both Ci and C2 in reflection one 
does not expect a very good agreement with theory, the problem is that the precise behavior of the diffuson near 
the surface is important. The diffusion approximation ^ 1 is not valid anymore, in particular if the indices of 
refraction match (then the intensity gradient near the surface is the largest). The calculation of Ci is probably easily 
extended, yet for C2 one also would need the Hikami box beyond first order in q£. 



D. Conductance Fluctuations: C3 

One could be tempted to think that the calculation of the C3 or UCF in the above Landauer approach is also 
straightforward. However, the calculation in the Landauer approach is quite cumbersome, as on short length scales 
(of the order of the mean free path) divergences show up when one treats the problem on a macroscopic level using 
diffusons. As an alternative, the Kubo approach is often used in mesoscopic electronic systems (Altshuler, 1985; Lee 
and Stone, 1985; Lee et ai, 1987). The results for the conductance obtained by either Kubo or Landauer formalism 
should be identical, see Fisher and Lee (1981) and Janfien (1991). Yet the Kubo approach cannot be applied directly 
to optical systems when absorption is to be included. Studies of C3 in the Landauer approach were undertaken by 
Kane et al. (1988), Berkovits and Feng (1994), and by the present authors (van Rossum, Nieuwenhuizen and Vlaming, 
1995) where it is shown explicitly that all divergences cancel. 



The C3 correlation function, defined in Eq. (344), involves incoming diffusons /Z"^*^, and outgoing ones C'^'^^. Due to 



the factorization of the external direction dependence, see Eqs. ( |349| ), ( |352| ), and (356), this dependence on external 
parameters cancels from C3. We can write 

Ct''{K,n)^Cs{^,n) ^^Y.F{q^,n,n). (369) 

where q_L is the {d — l)-dimensional transverse momentum. The function F is the main object to be determined 
in this section. We thus calculate it at fixed q± and with external diffusons for diffuse illumination iZin and £out, 



Eqs. (347) and ( |352| ). Depending on dimension we have 

F{q^, K, Cl) ^ F(0, K, tl) quasi 1-d, (370) 



qj- 



W I ^F{qj^,K,il) quasi 2-d, (371) 
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= W^ J ^F{qi_ , «, n) 3-d. (372) 

In contrast to Ci and C2, as a result of the q± integral, C3 will depend on the dimensionality of the system. 

Let us work out the long range diagrams of C3, roughly depicted in Fig. in more detail. For all diagrams there 
are two incoming advanced fields, which we momentarily term i and j, and two retarded ones, i* and j*. Consider 
an optical experiment with integrating spheres at incoming and outgoing sides. The incoming diffusons cannot have 
a momentum or frequency difference, and the pairing must be ii* and jj*. In Fig. ^(a) the incoming diffusons 
interfere somewhere in the slab. In a diagrammatic language the diffusons interchange a propagator so that the 
pairing changes into ij* and ji*. Propagation continues with these diffusons, which, due to the different pairings 
can have nonzero frequency difference and nonzero momentum. Only the outgoing diffusons without a momentum 
or frequency difference are dominant. Therefore, somewhere else in the slab a second interference occurs. Again 
exchanging an amplitude, the original pairing, ii* and jj*, is restored and the two diffusons propagate out, see 
Fig. p4|(a)i. Other contributions occur as well, yet the incoming and outgoing pairings are always ii* and jj* . In 
Fig. ^(a)ii the first incoming diffuson meets an outgoing diffuson and exchanges amplitudes. These internal diffusion 
lines meet at a second point where the original pairings are restored. Clearly in this process the intermediate paths 
are traversed in time reversed order. Due to time-reversal symmetry they give a contribution similar to Fig. p^(a)i. 
In Fig. |2j(b)i a diffuson breaks up such that one of its amplitudes makes a large detour, returns to the breaking point 
and recombines into an outgoing diffuson. The second incoming diffuson crosses the long path of the amplitude, and 
one of its amplitudes follows the same contour as the first diffuson. Finally, Fig. |2j(c) depicts the situation where 
only one internal diffuson occurs. Its endpoints must lie within a distance of a few mean free paths. Because of the 
local character this class does not show up in the final result; we need it, however, since it contains terms that cancel 
divergences from the other two classes. 



E. Calculation of the C3 correlation function 



The diagrams for the conductance fluctuations contain a loop; the two internal diffusons have a free momentum, 
over which one has to integrate. In Fig. ^(a) we denote this momentum q. Physically, one expects important 
contributions to the conductance fluctuations if the distance between the two interference vertices ranges from the 
mean free path to the sample size. Yet the q— integral diverges for large momentum, i.e. when the two interference 
processes are close to each other. The standard picture of diffuse transport with diffusons and interference described 



by Hikami vertices, that works so well for loop-less diagrams such as the C2 correlation function (section XV), and 



the third cumulant of the total transmission (section XVIl), becomes spoiled by these divergences 



As an example we calculate the diagram presented in Fig. |2^(a)i. This diagram was flrst depicted by Feng et al. 
(1988) and considered in detail by Berkovits and Feng (1994). These authors pointed out that a short distance 
divergency appears. They found a cubic divergency in three dimensions, and in general, a d-dimensional divergency 
in d dimensions. Let us see how it appears. For simplicity we consider a quasi one-dimensional system for which 
frequency differences and absorption are absent, therefore the decay rate vanishes, i.e. M = for all diffusons. From 
Fig. |2^(a)i one directly reads off its corresponding expression Fa.i 

Fa. = J J dzdz'Cl{z)H,{z)£l,{z,z')H4{z')Cl,,{z). (373) 

After performing some partial integrations and using the diffusion equation Eq.( |350D one finds 



Fa. = ^S{0) J dzCUz)Cl^,{z) 



+^ I dz/:i„t(z,^)[/:£(z)/:LW + 4>(^)>c?.t(^) 
+2/:;„(z)£;,t(z)£i„(z)/:out(z)] 

+Ahl f dz' f dzCUz,^')^'^iz)Kltiz'). (374) 



were the spatial derivative of C{z) was denoted The divergency is clearly present in the S{0), it comes about as 
follows: For the case of zero external momenta, the Hikami box yields iJ4(q, 0, — q, 0) = 2/14(7^, while the internal 
diffuson has the form £int(<z) = 127r/(^'^q^). Omitting the external hues, the diagram leads to 
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'^'^-i?|(q,0,-q,0)/:L(<j) = -?, / -f^q'^ = ^S{r^O). (375) 



We label the expressions according to the diagrams in Fig. 24, Fa, Fb and Fc and sum them. For the ID case we 
find 

Fa + Fb + F, = -^5{Q)L + ^. (376) 
15 15 



The second term in Eq.(376) is exactly equal to the well known result for the UCF in one dimension. (Note that the 
pre-factors of the diffusons and the Hikami boxes have canceled precisely. This is a manifestation of the universal 
character of conductance fluctuations.) But a singular part is annoyingly present. The term 5{Q) is a linear divergency. 
In the three-dimensional case one has to integrate over the transverse momenta yielding a cubic divergency. 

The cancellation of this divergence is involved and requires a number of special short range diagrams which resist 
a general classification or treatment (van Rossum, Nieuwenhuizen and Vlaming, 1995). Here we only use that finally 
the special short range diagrams cancel the divergences, yet they give no long range contribution. The nondivergent 
part of the diagrams of Fig. ^ is the only remainder. It reads 

F(gi,K,i7)=4/i2 j Jdzdz' Cint{z,z';M)Cintiz,z';M*) 

X [/:;f,(z)/:^2^,(z') + /:|„(z)/:^,,(z)/:|„(z')^out(^')] 

I Jazdz' [CUz,z';M)+CUz,z';M*)]^[CUz)Cout{z)] 
x-^[Ciniz')Coutiz')]. (377) 



in which again M = q± + k + ifl. The upper line of (377) corresponds to th e diagrams of Fig. |24|(a), whereas 



the lower line corresponds to the diagrams of Fig. ^(b). Finally, with Eq. (372) the value at vanishing transverse 
momentum gives the variance of the conductance in one dimension, whereas integration over the transverse momentum 
yields the correlation in two and three dimensions. 



Using Eqs. (|377| ) and ( 372 ) and inserting inserting the appropriate diffusons, various cases can be studied. Consider 



the case of fully transmitting surfaces; neglecting absorption and frequency differences one finds 

p^^^^ ^ 3 2 + 2q^^L^ - 2cosh2g^L + g^Lsinh2g^£ ^^^^^ 
2 qi_'^L'^ sinh^ q^L ' 

which decays for large momenta as q±~^, thus converging in three or less dimensions. The subsequent integration 
over the loop momentum finally yields 

{T^)c = ~ 0.133, quasi ID (379) 
15 

3 A A 

= ^C(3)- ~ 0.116-, quasi 2D (380) 
n-^ L L 

= i^-0.159A 3D, (381) 

in which ( is Riemann's zeta function. These values are for wide slabs. For cubic samples Lee et al. (1987) find in 
2D that (T^)c0.186 and in 3D a value of 0.296 (the value in ID is of course the same). One can also determine the 



frequency dependency of the correlation using (377), which is of importance as it determines the frequency range of the 
light needed to see the fiuctuations (van Rossum, Nieuwenhuizen and Vlaming, 1995). By inserting the appropriate 
diffusons one can also study the case of partial reflection at the surfaces of the sample. It is physically clear that the 
internal reflections lead to a less steep diffuse intensity in the sample as a function of the depth. The fluctuations are 
proportional to the space derivatives and thus reduce. We present the results in Fig. |2^ where we have plotted the 
correlation function for various values of the ratio between extrapolation length and sample thickness. One sees that 
the correlation is lower than without internal reflections. We note that neither the variance (the value at vanishing 
O), nor the form of the correlations are fully universal. 
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XVII. THIRD CUMULANT OF THE TOTAL TRANSMISSION 



So far we considered correlations between intensities defined as 

{TabTcd) 



{Tab) {Ted) 



= 1 + Ci + C2 + C3. (382) 



The Ci((7x, il), C2(9_L, and 6*3(17) correlation functions describe correlations between two intensities. By definition 
their value at zero frequency and zero momentum difference, i.e. the peak value, is the second cumulant of the 
distribution functions of the transmission quantities. An obvious question is how higher cumulants behave. The size 
of the fluctuations and the shape of the distribution relates to the "distance" from the localization transition. Far from 
localization, diffusion channels are almost uncorrelated and fluctuations are small (except the optical speckle pattern 
in the angular resolved transmission). The correlation between the channels increases if localization is approached. 
The relevant parameter is the inverse dimensionless conductance which can be interpreted as the chance that 
two channels interfere. Close to Anderson localization g approaches unity, and fluctuations increase. 



A. Cumulants of the probability distribution 

We first study the third cumulant of the total transmission. The total transmission is a constant superposed with 
fluctuations. To first order in the fluctuations have a Gaussian distribution (de Boer et ai, f992; Kogan, Kaveh, 
Baumgartner and Berkovits, 1993). The relative variance of this distribution, the top of C2, is proportional to g~^, it 
is thus a factor g larger than for the conductance fluctuations. This sensitivity of the total transmission to interference 
processes and its simple limiting behavior (as compared to the angular resolved transmission) make it an ideal quantity 
to study mesoscopic effects in its distribution. The third cumulant of the distribution was determined in very precise 
experiments by de Boer, van Rossum, van Albada, Nieuwenhuizen and Lagendijk (1994). In this section we calculate 
its theoretical value. Using the diagrammatic approach we relate the third cumulant normalized to the average total 
transmission, ((T^)), to the normalized second cumulant {{T^))- 

The moments of the probability distribution can be extracted from the distribution P{Ta) as 

{T^) = j ATaP{Ta)T^. (383) 

In a diagrammatic approach the fc-th moment can be represented by a diagram with k diffusons on both incoming and 
outgoing sides. The k = \ term is the average total transmission (Tq), as given by the Schwarzschild-Milne equation 
in Eq.( |355| ), it corresponds to a single diffuson and is thus independent of channel-to-channel correlations. The second 
moment can be decomposed in the first two cumulants: 

(T2) _{Ta)^ ^ {T^)cum ^^^^^^2^^^ ^33^^ 



(T,)2 (r,)2 (Ta) 

where the double brackets denote cumulants normalized to the average. Diagrammatically we depict the second 
moment in Fig. The decomposition in cumulants proves useful as each cumulant corresponds to a different 
number of interactions between the diffusons. The first term, Fig. |2^(a) contains no interference; it factorizes in the 
average transmission squared. The second term. Fig. ^(b), is the second cumulant {{T^))- It gives the variance of 
the fluctuations. This is just a special case of the C2 correlation function, ((T^)) — 6*2(0). 
Likewise, one can distinguish three different contributions to the third moment, 

^^1 + 3{{T^)) + {{t!)). (385) 

We have drawn the corresponding leading diagrams in Fig. J2^. The first term, Fig. p7|(a), again corresponds to the 
transmission without interference. The second term. Fig. ^(b), is the product of a single diffuson and a second 
cumulant diagram. From the figure it is clear that this decomposition can be done in three ways which determines 
the combinatorical prefactor of {{T^)) in Eq. ( |385|) . The third contribution stands for the third cumulant of the 
distribution and expresses the leading deviation from the Gaussian distribution. This is the term we are mainly 
interested in. It consists of two related diagrams: Fig. ^(c-|-d). The three intensities can interfere twice two by two, 
or the intensities can interact all three together, with a Hikami six-point vertex. Both contributions will prove to be 
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of the same order of magnitude. The strength of the effect can be easily estimated using the interpretation oi l/g 
as an interaction probabihty. By looking at the diagram, the third cumulant is proportional to the chance of two 
diffusons meeting twice, thus of the order \/g^. We have thus found the estimate 

{{TD) « (386) 

Finally, we note another that there is another contribution to both second and third cumulant. Because there is 
only a finite number of channels in an experiment, the intensity distribution will always have non-zero width. We 
have verified that this effect brings a negligible contribution to the measured second cumulant and the third cumulant 
(van Rossum, de Boer and Th. M. Nieuwenhuizen, 1995). 



B. The calculation of the third cumulant 



Two processes contribute to the third cumulant. One with two four point vertices, that we term ((T^))c and one 
with a six point vertex ((T^))^, where we have chosen the subscripts according to Fig. 



Interference via two four point vertices 

First consider the diagram in Fig. |2^(c). We have labeled to incoming diffusons with odd numbers, the outgoing 
ones with even numbers. Two incoming diffusons, and £3, meet at a position z, in a Hikami box they interfere 
into £2 and an internal diffuson £73*. £2 propagates out, whereas £7^* interferes again at z' with the incoming 
diffuson £5 into two outgoing ones, £4 and Cq. Apart from this process, also three other sequences of interference are 
possible. This means that the diffusons can also be permuted as: (£J, £3, £5, £2, £4, £6) (£3 , £5, £°, £4, £5, £2) 
(£5, £°, £3, £6, £2, £4)- We denote the sum over these permutations as X^per- ^^'^ diagrams can also be complex 
conjugated, there is also a combinatorial factor 2 for all diagrams. (Note that this is different from the C2 calculation; 
the second cumulant diagram is identical to its complex conjugate and thus there is no such factor.) The expression 
for the diagram of Fig. ^(c) is now 

((T3)), = (ra)-32VA Tdz r dz' Hi{z)Hi{z')Cl{z)C2{z)C:i{z) X 

per 

U{z')C5{z')C^{z')C)^i{z, z'). (387) 



In turns out that it is useful to rewrite the form of the Hikami boxes as introduced in Eq.(339) into an equivalent 
expression, using momentum conservation + (\a' + Qb + Qb' — 0. In real space the use of momentum conservation 
corresponds to partial integration. The Hikami box is again simplified using the fact that there are no transverse 
momentum terms, or qj^ terms, for the outgoing diffusons. Using the numbering in Fig. |2^, we obtain 

Hi{z) = ~hi[2d,,d,, + 29,,<9,3], Hi{z') - hi[2d,,d,, - 8^ + qi_l]. (388) 

Source terms, i.e., qf— terms of the incoming and outgoing diffusons were again neglected, but the source term of the 
diffuson between the vertices is important. As one sees with the diffusion equation ( |350D , it brings 

127r' 



i-dl + q4)^\^i{z, z') —5{z z'). (389) 



The contribution from the source term, i.e. H4,{z') oc — + qi_l, H4,{z) oc dzidz2 + SzjC^zg, is 



/fc2 



- / Az[dz,dz,+dz,dz,+ ... + dzMCl...C^. (390) 
Jo 



Although this corresponds to a local process (just one z coordinate is involved), it is of leading order. Together 
with the expression coming from H4{z'), proportional to dz^dz^, we find for the total contribution of the process in 
Fig.|^(c) 
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oi,2 i-L 



z) X 



per 

\z' C'^{z')C,{z')W)d,C}-\z,z') 





i A f 

x/:i(z)/:2(^)>c^(^)A(^)/:5(2)/:6(^), (39i) 



where C'{z) denotes the z derivative of C{z). Calculated for a plane wave, one finds 

{{t!))c^-^, (392) 
which is indeed proportional to g~^, as predicted. 

Contribution of the six-point vertex 

The other diagram contributing to the third cumulant is depicted in Fig. ^(d). The hexagon is again the Hikami 
six point vertex Hq. It can be thought of in the following way: the use of the Hikami box assumes that the outgoing 
legs scatter at least once before they propagate out or interfere again. This is a reasonable assumption for the outgoing 
diffusons, but for the internal diffuson £7^* it is also possible that coming from z it directly, i.e., without scattering, 
interferes again at z' . This process was not yet included in the calculation of the previous subsection but has to be 
studied separately. The unscattered intensity decays exponentially over one mean free path, therefore this process is 
only important if z and z' are within one mean free path. As the diagrams can also be complex conjugated, there is 
also a factor 2 for all diagrams. 

After a Fourier transformation in the z— direction the six-point vertex yields a contribution to the third cumulant 



3 _ hd^A 



^A f 

/ dz[d2^d,_^+d,^.,d,^^+dz3d,_^+d;:^dz^+dz,^dzg+d;,gdz^] 

a) Jo 



k^Ta 

xC,{z)C2{z)C3{z)C4z)C5{z)£e{z), (393) 



were again the notation of Eq. ( |342| ) is implied. We used the fact that all outgoing diffusons have zero transverse 
momentum. Therefore all qx^g^ -terms are absent. In case of an incoming plane wave we find a contribution to the 



third cumulant {{Tj^))d — — 4/(5g^). The contribution from the source term, i.e. Eq.(B90) of the previous subsection, 
exactly cancels the contribution from the six-point vertex. The cancellation seems plausible as one does not expect 
short distances properties to be important in the total process. Nevertheless, this cancellation depends on the precise 



form of the Hikami four-point vertex in Eq.(388 ). If we use other equivalent forms of the Hikami box the contributions 
of the single and double integral i n Eq .(391) shift with respect to each other and a full cancellation is not present. Of 
course, neither the result for Eq.(391) nor the final result for ((T^)) relies on this choice. One obtains for the third 



cumulant 



ShjAk^ 



^ / dz/:?(z)4(z)£^(z)x 



3 

per 
L 

dz' A(^')A(^')4(^')<9.£'"*(^, z'). (394) 



C. Influence of incoming beam profile 

Now that we know the leading interference processes, inserting the diffusons gives the final value of the third 
cumulant. We first consider the simple case of incoming plane waves. As there can be no transverse momentum 
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difference in the incoming amphtudes, all qJ_^ vanish. As a result all diffusons are linear functions of z. We find from 



Eq. (394) 



{{T!)) = ^ = y ((T^a))' plane wave; po » L. (395) 

In practice, however, one often deals with a Gaussian beam with limited spot size, influencing the cumulants in two 
ways. First, if the spot size decreases to values comparable to the sample thickness we have to convolute over a range 
of incoming momenta, just like we did when calculating the second cumulant. Second, the Gaussian profile brings 
an extra geometrical factor as will be shown below. The expression for the diagrams with diffusons with arbitrary 
momentum is calculated as follows. Because of momentum conservation and phase condition on the outgoing diffusons, 
the transverse momentum q_L7 of the diffuson connecting the two four-boxes must equal q_L5. The integration over 
the possible momenta results again in a Gaussian weight function exp{-_p^{q±l + q±^ + q±l)/8). The third cumulant 



is obtained by inserting the momentum dependent diffusons into Eq.(394). This gives 



{{T!)) = J d^q^id^q^s e-''"[«-?+'-3 + («.i+9.3)^]/8^3(^^^i^^^^i^ |q^l +q^3|L), (396) 



with 

F3{xi,X3,X5) 



y^r (xi +2:3)^X5 cosh(a::i +a:;3) (xi — 0:3)^X5 cosh(xi — X3) 

{xi+X3+X5)^{xi+X3 - 2:5)2 {xi - X3+X5)^{xi - X3- x^y 

(xi+a;3)a;5 sinh(2;i+X3) ^ {xi - X3)x5 sinh(xi — 2:3) 
(xi+2;3+a;5)(xi-|-X3 - 2:5) (xi - 2:3 + 2;5)(2;i - 2:3 - 2:5) 
(xi +2:3) cosh(a;i + 2:3 + 22:5) (2:1 +2:3) cosh(2;i -|-2;3 — 2x5) 



+ 



4(xi-|-X3-|-X5)2 4(xi+X3-X5)2 

(xi — X3) cosh(xi — X3 + 2X5) (xi — X3) cosh(xi — X3 — 2x5). 



4(xi-X3 + X5)2 4(^x1 - X3 - x^)'^ 

[x5 sinh(xi) sinh(x3) sinh^(x5)] . (397) 



We study again the behavior if the beam diameters are wide. In the limit of large beam diameter (po 3> L) one finds 
^3(0, 0,0) = if, this means for the third cumulant ((T^)) = 4^3(0, 0, 0)/3g^ or 

{{t!)) - ^ ((r2))2, (Gaussian profile; po » L), (398) 



which differs by a factor | from the plane wave limit Eq.(395). This is a geometrical effect, depending on the profile 



of the incoming beam. In a real space picture one understands this effect easiest: The correlation depends on the 
distance, it is strongest if the incoming beams are close together. Therefore it is not surprising to see the infiuence 
of their overlap. In next section we calculate this geometrical factor also for higher cumulants. For the experi men tal 
relevant case that the beam diameter is roughly equal to the thickness, it turns out that the behavior of Eq. (|39§| ) is 
found for a large range of beam diameters. The increase of the correlation for smaller beams turns out to be roughly 
the same for both the third cumulant and the second cumulant squared. 

In van Rossum, de Boer and Th. M. Nieuwcnhuizen (1995) theory was extensively compared with the experiment 
reported in (de Boer et al., 1994). Apart from the correction for the finite beam diameter, as discussed above, two 
other experimental corrections were included: Internal reflections and contributions from disconnected diagrams. All 



corrections to (398) turn out to be relatively small. By presenting the results as the ratio between the second cumulant 



squared and the third cumulant, errors in the sample thickness and the mean free path cancel. For the experimental 



data it was found that ((T^)) = (3.3 ± 0.6)((T2))2, in good agreement with Eq. ( |398| ). 

Stoytchev and Genack (1997) performed microwave scattering on systems with values of g approaching unity. In 
that case the plane wave prediction applies. They found {{T^)) — (2.38 ± 0.05)((r^))2, indeed in good agreement 
with Eq. ( ^95| ). 

The extension of the calculation to higher cumulants is straightforward. The n-th cumulant will contain [n — 1) 
Hikami four-point vertices. So the contribution is ((T")) oc g^~". But it is clear that the calculation becomes laborious 



at large n. We follow another approach to this problem in Sec. XVIII 
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D. Third cumulant of the conductance 



The reader will not be surprised to learn that the third cumulant of the conductance can also be calculated 
using a diagrammatic approach (van Rossum, Lerner, Altshuler and Nieuwenhuizen, 1997). This third cumulant has 
the particular behavior that it vanishes in one dimension(Macedo, 1994). We already saw that in the mesoscopic 
regime the conductance shows universal fluctuations (UCF). The conductance being a random variable showing 
large fluctuations, one should consider its full distribution. It was soon clear that the first higher cumulants of the 
conductance are proportional to (Altshuler, Kravtsov and Lerner, 1986) 

(5")cu,n (X n < go, {g) » 1. (399) 

Here go is the mean conductance at the scale i. In the metallic regime far from localization, where (17) ^ 1, the higher 
cumulants are probably small, and the distribution of the conductance is roughly Gaussian. However, for n ^ go the 
decrease in magnitude of cumulants as described by Eq. ( |399| ) is changed into a very rapid increase ((x eicp[gQ^n^]). 
This leads to log-normal tails of the distribution (Altshuler et ai, 1986). With increasing disorder, the log-normal 
tails become more important. The full conductance distribution at the threshold of localization {{g) ~ 1) is at 
present unknown, yet it is quite plausible that the whole distribution crosses over to a log-normal shape in the 
strongly localized regime (see Refs. (Shapiro, 1990) and (Altshuler, Kravtsov and Lerner, 1991) for a discussion). 
Indeed, it is well known that in the strongly localized regime in one dimension the conductance is given by the 
product of transmission amplitudes, yielding a log-normal distribution (Anderson, Thouless, Abrahams and Fisher, 
1980; Abrikosov, 1981; Mel'nikov, 1981; Economou and Soukoulis, 1981). 

Although higher order cumulants govern the tails of the distribution, they do not affect it near the center. A 
deviation from the Gaussian distribution near the center is revealed, first of all, in the lowest nontrivial cumulants. 
An important step in this direction was the calculation of the third cumulant of the distribution using random matrix 
theory by Macedo (1994). He found for the orthogonal (/3 = 1) and symplectic ensemble {§_= 4) that the third 
cumulant of the conductance is proportional to l/g"^, thus the leading term given by Eq. ( |399| ) vanishes (see for 
instance (Mehta, 1967) for the definitions of the ensembles). For the unitary ensemble {(3 = 2) even this sub-leading 
term vanishes, and the behavior is 1/g^ (Macedo, 1994). The physical reason behind this is not clear. (We recently 
learned that the same result in quasi one dimension was found Tartakovski (1996) by the scaling method described 
in (Tartakovski, 1995).) However, random matrix theory is only valid in quasi one dimensional systems. Therefore, 
it was not known whether such a cancellation also occurs in higher dimensions. If this is indeed the case, this might 
indicate an overlooked symmetry of the system. That this might be the case is also suggested by the fact that the 
leading order contribution to the third cumulant of the density of states vanishes in two dimensions (Altshuler et at, 
1986). The question we thus wish to consider is whether or not a cancellation occurs for the third cumulant of the 
conductance in two and three dimensions. 

We rewrite Eq.(399) into a relation for the relative cumulants: (5")c/(5)" oc The inverse powers of (g) on 



the right hand side can be interpreted as the number of Hikami four-boxes in the diagrams. Thus the third cumulant 
diagram contains four Hikami four-boxes. Indeed, it proves impossible to draw {g^)c diagrams with three boxes. 
Diagrams with more boxes are sub-leading as they are of higher power in 1/ (g) . The diagrams are represented in Fig. 
p8| . In the figure there are diagrams of the same order with one six-box and two four-boxes. They can be obtained 
by contracting one diffuson in the diagrams with four-boxes. In physical terms the six-box diagrams correspond to 
processes where after an interference process, the amplitudes do not combine into a diffuson. Instead, the amplitudes 
interact again directly without being scattered. 

The diagrams were evaluated using a Kubo approach, the advantage being that no divergences emerge; treatment 
of absorption is difficult though, as indicated in the discussion on C3. The evaluation of the diagrams is a painful 
but straightforward exercise. The diagrams have two or three free momentum loops, of which one is eliminated due 
to momentum conservation. In one dimension, where the diffusons are simple linear functions, the diagrams can be 
evaluated analytically. Numerical evaluation of the integrals in d = 2 and d — 3 yield 

{g^)c = (quasi ID), 

(c,3)^ = -0.0020(5) (quasi 2D, square), 

(c,3)^ = -f0.0076(.g)"i (3D, cubic). (400) 

Thus the leading contribution to the third cumulant in one dimension vanishes. This confirms the random matrix 
theory result (Macedo, 1994) diagrammatically. It is seen that there is no cancellation in higher dimensions. The 



results for rectangular samples are given in Fig. 29, where we multiplied the third cumulant by the average of the 
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dimensionless conductance. The third cumulant for wide slabs {L^ ^ L^, Ly ^ L^) is proportional to {LxLy/{g)). 
In the figure, due to the multiplication by (g), there is proportionality to (LxLy)"^ / L* for wide slabs. Note that the 
third cumulant passes through zero when going from 2D to 3D at a sample size of OAQL^ x L^y. L^- For very narrow 
slabs the correct quasi 2D limit is recovered. 

The random matrix result that the third cumulant vanishes in one dimensions is thus also found in the diagrammatic 
approach. In two and three dimensions, however, the cancellation is not present; the leading contribution to the third 
cumulant is negative in two dimensions and positive in three dimensions. The fact that the third cumulant changes 
sign is surprising. The third cumulant is also known as the skewness of a distribution. In analogy with the third 
cumulant of the total transmission, or if the distribution would be tending to log-normal, one would have expected 
a positive third cumulant of the conductance. Instead, we find that all possible values occur: negative, positive and 
zero. Experiments measuring the conductance distribution could enlighten this point. 



XVIII. FULL DISTRIBUTION FUNCTIONS 

In this section we calculate the full distribution of the total transmission and the angular resolved transmission; both 
can be mapped on the eigenvalue distribution of the transmission matrix. Before we discuss the distribution in the 
mesoscopic regime, let us first discuss the 'classical' situation, which is the limit of large g. The correlations between 
intensities are now absent. Here the angular resolved transmission, or speckle intensity, is distributed according to 
the Raleigh's law (Goodman, 1975), 

^'M = 7;^e-^-/<^->. (401) 

(Jab) 

Anyone who ever saw a laser speckle pattern on the wall will remember the wild pattern with both bright and dark 
spots. This is due to the Rayleigh's law. Rayleigh's law can be derived in the context of the ladder diagrams (Shapiro, 
1986). The n— th moment of the speckle intensities is made up of n amplitudes and n complex conjugated amplitudes. 
As stated before, there is no restriction on the pairing of two amplitudes into a diffuson for this type of measurement. 
Therefore the amplitudes {a^, ag, . . . , a*} have n\ possibilities to pair with the amplitudes {ai, 02, 03, . . . , a„}. 
Thus 

{T:,) = n\{Ta,r (402) 

which corresponds with Rayleigh's law. Yet, as we know from the example of a laser beam reflecting from a rough 
wall, multiple scattering is not essential for Rayleigh's law. Indeed, Rayleigh's original derivation goes as follows 
(Goodman, 1975): The field ip at a given position on the outgoing side is the sum of many fields, and therefore the 
real and imaginary part each have a independent Gaussian distribution. The intensity / is the amplitude squared and 
has thus an exponential distribution 

P ^ Q-m^i'f+C^^i'?]/'^'^'' ^ Q-^K^). (403) 

Thus reflection from a rough wall will also do. 

Kogan and Kaveh (1995) addressed the question of the speckle statistics of (optically) very thin samples. In optically 
very thin samples a considerable part of the light does not scatter at all, but is transmitted coherently. In the limit 
of zero thickness, i.e. only coherent light, the distribution becomes a delta-function. In the intermediate regime a 
simple counting argument yields the distribution, which interpolates between the delta distribution and Rayleigh's 
law. In the opposite case of strong scattering, the mesoscopic regime, interferences modify the speckle distribution. 
The leading correction was derived by Shnerb and Kaveh (1991). Genack and Garcia observed a deviation from 
Rayleigh's law at large intensities (Genack and Garcia, 1993; Garcia and Genack, 1989). A crossover to stretched 
exponential behavior was derived by Kogan et al. (1993). 

For the total transmission and the conductance the distributions are similar in the large g limit. For both quantities 



the outgoing diffusons are frequency and momentum independent as we explained in Sec. XV. As there is no interaction 
between the diffusons in this limit, only one type of pairing is possible for both quantities. This yields in principle a 
delta-distribution if one probes an infinite number of channels. In practice of course only a limited number of channels 
is probed, and the law of large numbers predicts a narrow Gaussian distribution. 
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A. Eigenvalues of the transmission matrix 



In mesoscopic systems the observables are random quantities and are therefore not always characterized by their 
mean values, but their entire distribution function is of interest. This is particularly prominent in the distribution 
of eigenvalues of the transmission matrix. Assuming that all eigenvalues contribute equally to the conductance, one 
expects a Gaussian distribution. As there are N eigenmodes, the eigenvalue distribution should be peaked around 
g/N = £/L. But this picture proves wrong, as was first put forward by Dorokhov (1984) and later also by Imry (1986). 
Not the eigenvalues but the inverse localization lengths 1/^n are uniformly distributed, see Pendry, MacKinnon and 
Pretre (1990) and Stone et al. (1991). As a result the eigenvalues have a "bimodal" distribution peaked around and 
1. The eigenvalues T„ of the transmission matrix tH can be expressed as 

T„ = cosh-2(L/^„), (404) 

implying that 

(Tr(ttt).-) ^ Ti) = g ,^^^j—T^- (405) 

This distribution is plotted in Figjs^. The first derivations of the distribution used random matrix theory and were 
therefore valid in quasi-lD only (Mello, Pereyra and Kumar, 1988). However, Nazarov (1994) showed that it is not 
only true in quasi-lD, but under very general conditions. Note that the normalization of the distribution is ill-defined, 
the distribution should be understood in the sense that all its moments are defined. This problem can be avoided by 
adjusting the lower boundary of the integral instead of 0. The minimal value correspon ds to the decay of unscattered 
intensity, therefore T„ > cosh~^(i/^) « exp(— 2i/^), normalizing the distribution Eq. ([405| ). Although this cutoff is 
important for the normalization of the distribution, its influence on the momenta is very small and we refrain from 
it. The distribution Eq. ( [405| ) is valid only in the regime where g is not too small. Loop effects near the Anderson 
transition, where 5 ~ 1, will change it, as we will discuss below. 

The concept of open and closed channels explains the physical meaning of this curious distribution function. An 
eigenmode of the transmission matrix is, according to this distribution, either essentially blocked or, with a much 
smaller probability, essentially conducting. This picture was confirmed in various computer simulations (Pendry et al., 
1990; Pendry, MacKinnon and Roberts, 1992; Oakeshott and MacKinnon, 1994). The eigenvalue distribution also 
provides us with a nice picture explaining why the correlations in mesoscopic systems are so large. If all the channels 
are equally conducting, fluctuations of the channels average out by the law of large numbers. Yet if only a few channels 
are conducting, fluctuations in one channel will be clearly seen. Using these concepts the probability distribution of 
the total transmission and deviations from Rayleigh's law will be explained below. 

We should point out that for the UCF the same argument goes wrong. The variance in the conductance Nl/L is 
caused by the finiteness of the number of channels N . In the picture of either closed or open channels the number 
of open channels A'eff equal roughly Ni/L <C N . The conductance would thus fluctuate according the law of large 



numbers with a variance: (^„ T^) = ^g. However, this is incorrect as we know that (^„ T^) ~ 0(1), see Sec. XVI. D 



3 

For the UCF the interaction between the eigenmodes is essential. 

In practice it is difficult, if not impossible, to measure the eigenmodes and eigenvalues directly. As they are the 
eigenmodes of the very large random transmission matrix, they have a very complex structure. Nevertheless, the 
eigenvalue distribution has observable consequences. It was shown by Beenakker and Biittiker (1992) that the shot 
noise of electronic conductors universally reduces by a factor 3 because of this distribution. Shot noise only occurs 
in electronic systems and is thus not of much interest for us. Instead, we will show that the total transmission and 
speckle intensity distribution function are related to the eigenvalue distribution function. 



B. Distribution of total transmission 



We first calculate the probability distribution of the total transmission. From the distribution without interference 
effects we saw that it has a simpler distribution (in the sense of cumulants) than the speckle distribution. We neglect 
corrections as absorption, skin layers, and disconnected diagrams. Consider again an incoming plane wave in direction 
/ia (/i denotes again the cosine with respect to the z axis). The wave is transmitted into outgoing channel b with 
transmission amplitude tab and transmission probability Tab = l^afcP- The speckle and total transmission factorize as 
{To) = Eaff and {Tab) = Caebff- 
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We consider the j*^ cumulant of T^. In a diagrammatic approach this object has j transmission amphtudes tab and 
an equal number of Hermiti an con jugates tl^ = t*^. The exphcit calculation of the second cumulant C2 (Sec. |XV| ), 



and third cumulant (section |XVI]| ) showed that the leading diagrams are connected and have no loops. 

Let us fix the external diffusons in the term ^objibia^aba ' ' '^ab 4 o- Contributions to the sum over bi only come 
from diagrams with outgoing diffusons that have no transverse momentum. These are the diagrams where the lines 
with equal bi pair into diffusons. We indicate this pairing of the outgoing diffusons with brackets 

itabAjitabAa)---itabA,a)- (406) 

The outgoing diffusons are now fixed. At the incoming side it is convenient to keep for instance the tatS fixed, while 
the ijj^'s are permuted among them. This leaves j! possibilities for the incoming side in total, but of these only 
(j — 1)! permutations correspond to connected diagrams (these diagrams contain j — 1 Hikami boxes. The C2 diagram 
contained one box, the third cumulant diagram two boxes) . Next we factor out the incoming and outgoing diffusons 
and group the remainder of the diagrams into a skeleton K. Making use of Eq. (p4^ we obtain: 



(7^i)con=ei(j-l)! J dridr'i • • • dr.dr; A„(ri)/:out(rl) 

• • • -Cin (r, )Cont {r'j)K{Yi , r'l , • • • , r„ r^. ). (407) 

The integral just describes 

(Tr(ttt)^) ^ Tr ^ KbAaMb.'-'^a^bA^a.) 

bi ,a2,b2,---,aj ,bj 

= X! (^aibl4ia2*a262 " ■ ^aj&j^Ijai)- (408) 

ai ,61 , ■ ■■,aj ,bj 

There is only one way to attach incoming and outgoing diffusons to K. The sums over the indices lead exactly to the 



total flux diffusons in Eq. ( 407 ) . One finds the following important relation between the moments of the eigenvalue 



distribution and the connected total transmission diagrams 

(T^')con = (j-l)!4(Tr(ttt)^). (409) 

Normalizing with respect to the average, we introduce Sa = Ta/{Ta). The generating function of the connected 
diagrams is easily calculated 

$con(a;) = 2^ (s^Jeon 

i=i 

= 5 log' (v/l + x/<?+ ■ (410) 

Since the cumulants are solely given by connected diagrams, the distribution of follows as (Nieuwenhuizen and van 
Rossum, 1995) 



r°° dx 

P{sa)^ — exp[a;s„ - $con(a;)] . (411) 

J_ioo 27^^ 



Let us examine some properties of the distribution. For limiting values of Sa we use a saddle point analysis. The 
saddle point is found by the condition ^ [xSa ~ $con(a;)] = 0. Thus we find 



log {^^X^xjg + v/aTs) 
\Jxlg{\ + xlg) 



(412) 



The r.h.s. of Eq.( [412| ) diverges if x approaches —g (from above) and decreases monotonically for larger x. Thus for 
large Sq 3> 1, we find the saddle point near x = —g. By inserting the saddle point one finds a simple exponential tail 

P{Sa) « exp(-5Sa + .g — ), Sa > 1. (413) 



67 



The saddle point also dominates the shape for small Sq <C 1 (and large g). One finds essentially a log- normal 
growth: 



exp 



log log- 



(414) 



In Fig. ^ the full curves depict the distribution Eq. ([41l| ) for some values of g. At moderate g we clearly see 



the deviation from a Gaussian. In interesting, recent experiments on quasi ID microwave scattering, Stoytchev and 
Genack (1997) measured the total transmission distribution for conductances as low as g ~ 3. Although one would 
expect important loop contributions for such low values of g, excellent agreement with the above prediction is found. 
To obtain the fit, the value of g had to be extracted from the variance of the distribution rather than from the measured 
conductance. This adjustment is probably necessary because of the strong absorption occurring in waveguides. The 
effect of absorption was calculated by Brouwer (1998) for a waveguide. For strong absorption the total transmission 
has a log-normal distribution, while the ratio of ((r^))/((r^))^ becomes 3.0 instead of 2.4. 



C. Influence of beam profile 



Above we have considered the case of an incoming plane wave. Again in optical systems a Gaussian intensity 
profile is more realistic. For the third cumulant we saw a nontrivial dependence on the incoming beam profile, 
suggesting that also higher cumulants are sensitive to this effect. For perpendicular incidence the incoming amplitude 
is ■0in(r) = (?!)((7a)V'in('')' where ^/'(^^(r) is the plane wave of Eq. (147), and where 



(/)(ga) = V2^poexp{-^plql). 



(415) 



We consider the limit where the beam is much broader than the sample thickness {po 3> L) but still much smaller 
than the transverse size of the slab (po ^ W). With a smaller beam diameter, the incoming transverse momenta, 
which are of order l/po: a-re of the order of l/L. The diffusons will then become the well-known cosh functions, 
Sec. XVI. A. Here the momentum dependence of the diffusons can be neglected; apart from a geometrical factor, the 
diffusons are identical to the plane wave case. Due to the integration over the center of gravity, each diagram involves 
a factor ^(Ssg.Sg'- In the j*'' order term there occurs a factor 



= / d'p\cl^{p)f^. 



For a plane wave we have |(/'(/3)| = \/A, and Fj = . For our Gaussian beam we obtain 

' J V 2 



(416) 



(417) 



It is thus convenient to identify Aq ~ ^t^Pq with the effective area of a Gaussian beam. (This definition is different 
from the one used in previous sections, were Aq — tt^q). As compared to the plane wave case, the j^^ order term is 
smaller by a factor 1/j for a Gaussian profile. This implies for the generating function of the connected diagrams 



$con(a;) 




(418) 



For small Sa (and large g), there is again a log-normal saddle point. For large Sa, the dominant shape of the decay is 
given by the singularity at a; = —g and again yields P{sa) ~ exp(— gSa). The shape of the distribution is quite similar 
to the plane wave case. Expanding the generating function we recover the results for the second and third cumulant 
obtained in Sees. XV and XVII. 
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D. Speckle intensity distribution 



We apply the same method to obtain the distribution of the angular transmission coefficient. The angular and 
total transmission distributions can be related to each other in a fairly simple way, because the interference processes 
are dominated by the same type of diagrams: the loopless connected diagrams. The different distribution for total 
transmission and speckle is the consequence of a counting argument only. In the plane wave situation the average 
angular transmission reads {Tab) = £a£bg- Let us count the number of connected loopless diagrams that contribute to 
'^ab ~ ^ab^la^ab ' ' ' ^ab^ba' uot Only masslcss outgoiug diffusous, but all pairings into outgoing diffusons contribute. 
This yields an extra combinatorial factor jl in the j^^ moment: 

racon=j!(j-l)!4e^,(Tr(tti)-'). (419) 

For the normalized angular transmission coefficient Sab — Tab / {Tab) , we define the generating function of the connected 
diagrams as (Nieuwenhuizen and van Rossum, 1995) 



with <I>con given by Eq. (410) for plane wave incidence and by Eq. (418) for a broad Gaussian beam, respectively. In 
contrast to the total transmission distribution, the cumulants are not only given by the connected diagrams. Kogan 
et al. (1993) showed that the extra summation of the disconnected diagrams corresponds to an additional integral, 
which finally yields a Bessel function: 

da; 



P{sab)= — i^o(2V=s:^)exp(-$con(a;)). (421) 

J-ioo 

The integral can be evaluated numerically. The speckle intensity distribution is shown in Fig. ^ for an incoming 
plane wave. 

For large g and moderate Sah, one has <i>con(a^) ~ x and we recover Rayleigh's law: P{sab) — exp(— Sq^)- The leading 
correction is found by expanding in l/g 



P{Sab) 



l + :^(sab-4s,b + 2) 



(422) 



A similar equation was derived previously by Shnerb and Kaveh (1991). 

For obtaining the tail of the distribution, one can again apply steepest descent, which yields 

P{sab) - exp(-2V5i^ ). (423) 

This stretched exponential tail of the distribution was also seen in experiments. Already in 1989 Garcia and Genack 
observed it in microwave experiments. But unfortunately their dynamic range was rather small. The fit of stretched 
exponential is therefore rather unprecise. The maximum intensity was 5 times the average where one does not 



expect to see the tail behavior, rather one observes a cross-over behavior. Using Eq. (418) one easily derives the 
speckle distribution due to an incoming beam with Gaussian profile; it leads to a different distribution with the same 
asymptotic behavior. 

E. Joint distribution 

It turns out that the speckle distribution and the total transmission distribution of a certain incoming direction are 
related. This is not surprising as a large total transmission for a given incoming direction will also be reflected in the 
individual speckles. The moments of the joint distribution are 

{4bsl)con ^{k + l-l)\ k\ (Tr(ttt)fe+'). (424) 

The combinatorial factor fc! is for the speckle pattern, the factor {k + I — 1)\ is the number of possible pairings of 
incoming beam. Defining aab as Sab = o'abSa, one obtains {c\},s"'') — fc! (to — 1)! (T™) with m ~ k + I. Thus 
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P(<^ab,Sa) = eyi]i[—(Jab)P{,Sa)- This shows that Oab and Sa are independent variables. Nevertheless, Sab and Sa are 
dependent. Their joint distribution is 

P{sab, Sa) = / 7r~ exp(a;Sa - <^con{x)). (425) 

Sa J 2m 

Integration over Sa or Sab yields again the distribution P{sab), or P(sa), respectively. This is somewhat surprising. 



Indeed, as the right hand side of Eq. (425) contains only connected diagrams, the same holds thus for the left hand 
side. This seems to contradict our earlier remark that P{sab) contained both connected and disconnected diagrams. 
However, in deriving Eq. ( |425| ) we have tacitly performed analytic continuation for negative I values. By performing 
the integration over Sa we find the correct expression for P{Sab)- 

In conclusion we derived the full distribution function of both total transmission and angular transmission coefficient 
by mapping it to the distribution of eigenvalues of Taf,. But at least for the lowest cumulants one can use this mapping 
the other way around. First, experiments on the distribution functions thus confirm the eigenvalue distribution. 
Secondly, from our detailed evaluation of the low order diagrams, i.e. the second and third cumulant, insight in the 
eigenvalue distribution can be obtained. It is clear from those calculations that also the eigenvalue distribution is 
based on loopless connected diagrams. 

The effect of absorption has been neglected here. Our explicit calculation of second and third cumulants and 
the experimental data show that the effects of absorption and skin layers can be incorporated to a large extent by 
considering the normalized cumulants, such as (T^)c/(T'a )c- Indeed, also in the presence of strong absorption it was 
seen that these dimensionless ratios are similar to the values calculated without absorption (Stoytchev and Genack, 
1997). A theoretical analysis of absorption effects was carried out by Brouwer (1998). 

The distribution of the conductance cannot be calculated using the above technique. The correlations in the 
conductance correspond to correlations in the eigenvalues. These correlations were not included, as we only used the 
distribution of a single eigenvalue, and not the joint distribution of multiple eigenvalues. 
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APPENDIX A: LOOP INTEGRALS 

The same loop momentum integrals recur often in the diagram calculations: 

h,i^ /(^G'=(p)G*'(p), (Al) 
where G{p) = [p^ — fc^ — nt\~^. For the simplest integral, 7i^i we find 

= / |^o(p)«-(p) = £ ^ v^/)v^^'y <*^' 

where /x^ = —k'^ — nt. The sum of the residues yields l/[47r(^ + 71)], and without absorption the optical theorem gives 

nlmt 1 , . „, 

t^ + -p=^ = j. (A3) 



Therefore 



= (A4) 
The general Ik,i integral can be found from the 7i,i integral using 

= 2^i;i^''' ' = 2^^^'=''- ^^^^ 

For example, 7i,2 = --^2,1 = -W^/iSnk) and 72,2 = £^/{8nk'^). We shall also use 

/ ^iP-^fG\p)G*\p) = lkVlk,i, (A6) 
where the factor 1/3 comes from the angular integral. If absorption is present n+Ji = [£{1 — k^£^/3)]~^, thus 

' 47r(/x + 7t) 47r V 3 y ^ ' 

In general there is a prefactor 

7„,„(x(1-k2^V3)"+™-i, (A8) 

where we have assumed that k£ 1. However, this prefactor needs only to be included for terms not proportional to 
q or tt, as we work to first order of {M£). 
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TABLE I. Numerical values of several fundamental quantities in radiative transfer theory for typical values of the index ratio 
m. 



m 


To 


n(l) 


7(1,1) 






O.Uo 


Zi. ( 


zi.o 


U.ioo 


3/2 


2.50 


10.8 


10.6 


0.269 


4/3 


1.69 


8.34 


7.94 


0.343 


1 


0.710446 


5.03648 


4.22768 


1/2 


3/4 


0.815 


5.39 


4.63 


0.479 


2/3 


0.881 


5.60 


4.85 


0.465 


1/2 


1.09 


6.25 


5.55 


0.427 



TABLE II. Anisotropic versus isotropic scattering. Toltr is the thickness of a skin layer; ri(l) and 7(1, 1), respectively, yield 
the transmitted and reflected intensities in the normal direction; B(fS) is the peak value of the enhancement factor at the top 
of the backscatter cone; Ttr-AQ = k\itr^O is the dimensionless width of this backscatter cone. The third row gives the relative 
difference in the two cases. 





TO Tl(l) 7(1,1) B(0) Ttr^Q 


Isotropic(rir = 1) 
Very anisotropic(Tfr » 1) 


0.710446 5.036475 4.227681 1.881732 1/2 
0.718211 5.138580 4.889703 2 0.555543 


Difference (%) 


1.1 2.0 15.7 6.3 11.1 



EDITOR: table 3 can be set narrower, I couldn't figure out how to make Latex do so. 



Object 


Small radius behavior 


Large radius behavior 


Absorbing sphere 
Transparent sphere 
Reflecting sphere 

Absorbing cylinder 
Transparent cylinder 
Reflecting cylinder 


~ 4 

Jieff «a£(f)3^ e"^, owl.20 
(1 + f) 

p ^ 3Re 


Q^R{1-^) , TO « 0.710446 
P«-J?3(l-2mii^^ roi« 0.4675 

(1 + 1^) 
i?eff « i? (1 - 2^) , TO « 0.710446 

Px-R^{l-^^), roiw 0.2137 



TABLE III. Capacitance Q and polarizabiUty P of small and large spheres and cyHnders immersed in an opaque medium. In 
column for large radius the diffusion approximation is given by the term outside the brackets. For narrow and wide cylinders 
the effective radius and the polarizability are presented. 



75 



FIG. 1. Schematic representation of multiple scattering in a slab. A plane wave, impinging at an angle 6 with respect to the z 
axis, is refracted to an angle 6' inside the multiple scattering medium, while it is partly specularly reflected. After the diffusive 
transport in the bulk, a wave arriving under angle 9'i, at one of the boundaries, partially leaves the medium under angle 6b, and 
is partially reflected internally. 
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FIG. 2. Solution of the Schwarzschild Milne equation as taken from table 17 in Van der Hulst (1980) for a slab of thickness 
L = 4£ with no internal reflection, a — 1. The precise solution near the incoming plane strongly depends on /i = cos 9', where 
9' is the angle of the incoming beam. Apart from a multiplication factor, the bulk behavior is the same; all lines extrapolate 
to L + 00 (thin dashed line). 
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FIG. 3. The dressed Green's function. 
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^^^m = + — Z — + — Z — E— + -Z — Z — E— + • • • 

FIG. 4. Average Green's function (heavy line): Thin lines, the bare Green's function; dashed lines indicate that scattering 
takes place from the same potential. The self-energy E only contains irreducible diagrams. It generates the amplitude Green's 
function according to the lowest line. 
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FIG. 5. Various intensity diagrams. (A) The ladder diagrams. Upper line: A typical ladder diagram before averaging. The 

crosses are scattering potentials, the Green's functions are bare. Lower line: Averaging over the possible scattering diagrams 
leads to the ladder diagrams, here a diagram with three common scatterers, where the circles represent t matrices and the 
Green's functions are dressed. The ladder sum contains diagrams with an arbitrary number of common scatterers. (B) 
Some scattering diagrams that are not elements of the ladder diagrams. The upper diagram is a maximally crossed diagram, 
responsible for the enhanced backscatter cone. The lower diagram is of second-order in the density. 
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FIG. 6. Angular resolved reflection coefficient for perpendicular incidence on a semi-infinite medium, for different values of the 
index ratio m. For m < 1 no radiation can exit the random medium at an angle exceeding the Brewster angle. 




FIG. 7. Angular resolved transmission coefficient for perpendicular incidence on a thick slab. The ratio of refractive indices, 
m, has been indicated. 
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FIG. 8. Handle diagram with three scatterers. The set of handle diagrams with an arbitrary number of scatterers generates 
the minimal set of self-energy diagrams necessary for flux conservation with the backscatter cone taken into account. 
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FIG. 9. Maximally crossed diagrams 
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Q 

FIG. 10. Backscatter cone of a semi-infinite medium at normal incidence for several values of m. The diffusive background was 
substracted. 



FIG. 11. At largo angles the backscatter cone is mainly determined by second-order scattering. The term with internal reflection 
has a larger optical path and therefore decays faster. 
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FIG. 12. Co mpa ris on o f exact values and diffusion approximation for ti(1), 7(1, 1) and Stq. The solid line is 4/T, the dashed 
line is Eqs. (M2) + (304). Even for m close to unity the asymptotic results are quite good. 




FIG. 13. Plot of the size factor Q/Qdiff of the capacitance of an absorbing sphere, against the ratio R/i. Full line: outcome of 
the numerical analysis; Dashed lines: small-radius and large-radius behaviors of Table 3. 



80 




2 4 6 8 
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FIG. 14. Plot of the size factor pjp'^^^ of the polarizabihty of a transparent sphere (lower curves) and a reflecting sphere 
(upper curves). Pull lines: outcome of the numerical analysis; Dashed lines: small-radius and large-radius behaviors of Table 3. 
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PIG. 15. Hikami four-point vertex. It describes the exchange of amplitudes of two incoming diffusons 1 and 3 into two outgoing 
diffusons 2 and 4. The dots linked with the dashed line denote an extra scatterer on which both amplitudes scatter. The solid 
lines are dressed amplitude propagators. 
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3 4 

PIG. 16. Six-point vertex Hq, describing the interaction of six diffusons. We did not draw possible rotations of the three right 
diagrams. In total there are sixteen diagrams. 
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FIG. 17. Four-point vertex, beyond the second-order Born approximation. Instead of three, there axe now eight diagrams 
to be calculated. The resulting expression, however, is apart from a renormalization of the mean free path the same as in 
second-order Born approximation. 





FIG. 18. First-order correction to the conductivity in the Hikami formalism. At the right we have written out the vertices in 
the second-order Born approximation. On the second line the first diagram in the r.h.s. of the upper line is drawn in a different 
way, showing it is a loop effect. 
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FIG. 19. Schematic picture of the different correlation functions present in the transmission of two intensities. The arrows 
distinguish between advanced and retarded propagators. Propagators with equal transverse momentum or frequency have 
the same style of line (dashed or solid). The black box contains in principle all contributions to (TabTcd)- The main part 
just factorizes, but correlations are present. The Ci correlation involves a simple interchange of amplitudes. The C2 and C3 
correlation involve a Hikami box (circled H); they are a perturbative effect for weak scattering. 




FIG. 20. A contribution of order to the C2 correlation. The boxes are Hikami boxes, the parallel lines are the diffusons. 
By following the amplitudes, one can check that the incoming pairings ij* and ji* change on the outgoing side into pairings ii* 
and jj* . This diagram is thus not contributing to the C3, which is also of order 1/g but requires diagrams without amplitude 
exchange. The exchange in pairing of incoming and outgoing amplitudes means that the diagram is a higher order C2 diagram. 
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FIG. 21. Ci angular correlation function plotted against the scaled perpendicular momentum difference: Solid line, small zq, 
no absorption; short dashed line, absorption {k — 2/1/), small zq; long dashed line, large skin layers (2:0 ~ L/3), no absorption. 



FIG. 22. Diagram of the long-range C2 correlation function. The shaded box denotes the Hikami four-point vertex depicted in 
Figs, na and O. 
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FIG. 23. F2, the normalized form of the angular co rrela tion function C2, plotted against the scaled perpendicular momentum 
difference: Solid line, small zq, no absorption Eq. (363); short dashed line, absorption (k = 2/L), small zq; long dashed line, 
large skin layers (20 = L/3), no absorption. 
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FIG. 24. Leading contributions to the conductance fluctuations apart from short-distance processes. The incoming diffusons 
from the left interfere twice before they go out on the right. The close parallel lines correspond to diffusons; the shaded boxes 
are Hikami vertices; q denotes the free momentum which is to be integrated over. 
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FIG. 25. The C3 frequency correlation in 3D as a function of frequency difference; different internal reflection values: upper 
curve, zo = 0; middle curve, zo = L/10; lower curve zo = L/5; no absorption. 
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FIG. 26. The two contributions to the second moment of the total transmission, (a) Independent transmission channels. This 
process is of order unity and is almost completely reducible to the mean value squared, (b) Two interfering channels; this is 
the second cumulant or C2 diagram of order 1/g. The close parallel lines are diffusons. 
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FIG. 27. The three contributions to the third moment of the total transmission, (a) independent transmission channels; it is 
of order 1. (b) correlation which can be decomposed into the second cumulant; this is of order g~^. (c) and (d) contributions 
to the third cumulant, 0{g~^). 
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FIG. 28. Set of diagrams for the third cumulant of the conductance. The wavy hnes denote diffusons (in the Kubo approach). 

Third cumulant of the conductance 

3D geometry 



<g >x<g> 



0.03 



0.02 



0.01 



0.00 




0.0 1.0 2.0 3.0 

width ratio: LJL^ 



FIG. 29. Third cumulant multiphcd by the average conductance in a 3D sample, plotted against the transverse size Lx- The 
geometry is x x Lz for the solid line. Note that the sign changes in going from a quasi- 2D to a 3D sample. The dashed 
line corresponds to the geometry Lx x x Lz- 
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FIG. 30. The bimodal eigenvalue distribution. The average value occurs only with a small probability, Instead the eigenvalues 
are almost all either or 1. 




FIG. 31. Distribution of the total transmission vs normalized intensity Sa, theory and microwave experiments for three different 
sample sizes (a,b,c) (Stoychev and Genack, 1997). For curve (c), the conductance is the lowest at 3.06 (absorption corrected). 
Reproduction with kind permission of the authors. 
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